Source code for microstructpy.geometry.n_sphere
# --------------------------------------------------------------------------- #
# #
# Import Modules #
# #
# --------------------------------------------------------------------------- #
from __future__ import division
import numpy as np
__author__ = 'Kenneth (Kip) Hart'
# --------------------------------------------------------------------------- #
# #
# NSphere Class #
# #
# --------------------------------------------------------------------------- #
[docs]class NSphere(object):
"""An N-dimensional sphere.
This class represents a generic, n-dimensional sphere. It is defined by
a center point and size parameter, which can be either radius or diameter.
If multiple size or position keywords are given, there is no guarantee
whhich keywords are used to create the geometry.
Args:
r (float): *(optional)* The radius of the n-sphere.
Defaults to 1.
center (list): *(optional)* The coordinates of the center.
Defaults to [].
radius : Alias for ``r``.
d : Alias for ``2*r```.
diameter : Alias for ``2*r``.
size : Alias for ``2*r``.
position : Alias for ``center``.
"""
# ----------------------------------------------------------------------- #
# Constructor #
# ----------------------------------------------------------------------- #
def __init__(self, **kwargs):
if 'r' in kwargs:
self.r = kwargs['r']
elif 'radius' in kwargs:
self.r = kwargs['radius']
elif 'd' in kwargs:
self.r = 0.5 * kwargs['d']
elif 'diameter' in kwargs:
self.r = 0.5 * kwargs['diameter']
elif 'size' in kwargs:
self.r = 0.5 * kwargs['size']
else:
self.r = 1
if 'center' in kwargs:
self.center = kwargs['center']
elif 'position' in kwargs:
self.center = kwargs['position']
else:
self.center = []
[docs] @classmethod
def best_fit(cls, points):
"""Find n-sphere of best fit for set of points.
This function takes a list of points and computes an n-sphere of
best fit, in an algebraic sense. This method was developed using the
a published writeup, which was extended from 2D to ND. [#bullock]_
Args:
points (list, numpy.ndarray): List of points to fit.
Returns:
NSphere: An instance of the class that fits the points.
.. [#bullock] Circle fitting writup by Randy Bullock,
https://dtcenter.org/met/users/docs/write_ups/circle_fit.pdf
""" # NOQA: E501
# convert points to numpy array
pts = np.array(points)
n_pts, n_dim = pts.shape
if n_pts <= n_dim:
mid = pts.mean(axis=0)
rel_pos = pts - mid
dist = np.linalg.norm(rel_pos, axis=1).mean()
return cls(center=mid, radius=dist)
# translate points to average position
bcenter = pts.mean(axis=0)
pts -= bcenter
# Assemble matrix and vector of sums
mat = np.zeros((n_dim, n_dim))
vec = np.zeros(n_dim)
for i in range(n_dim):
for j in range(n_dim):
mat[i, j] = np.sum(pts[:, i] * pts[:, j])
vec[i] += np.sum(pts[:, i] * pts[:, j] * pts[:, j])
vec *= 0.5
# Solve linear system for the center
try:
cen_b = np.linalg.solve(mat, vec)
except np.linalg.linalg.LinAlgError:
cen_b = pts.mean(axis=0)
cen = cen_b + bcenter
# Calculate the radius
alpha = np.sum(cen_b * cen_b) + np.trace(mat) / n_pts
R = np.sqrt(alpha)
# Create the instance
return cls(center=cen, radius=R)
# ----------------------------------------------------------------------- #
# String and Representation Functions #
# ----------------------------------------------------------------------- #
def __str__(self):
str_str = 'Radius: ' + str(self.r) + '\n'
str_str += 'Center: ' + str(tuple(self.center))
return str_str
def __repr__(self):
repr_str = 'NSphere('
repr_str += 'r=' + repr(self.r) + ', '
repr_str += 'center=' + repr(tuple(self.center)) + ')'
return repr_str
# ----------------------------------------------------------------------- #
# Equality #
# ----------------------------------------------------------------------- #
def __eq__(self, nsphere):
if not hasattr(nsphere, 'r'):
return False
if not np.isclose(self.r, nsphere.r):
return False
if not hasattr(nsphere, 'center'):
return False
c1 = np.array(self.center)
c2 = np.array(nsphere.center)
if c1.shape != c2.shape:
return False
dx = np.array(self.center) - np.array(nsphere.center)
if not np.all(np.isclose(dx, 0)):
return False
return True
def __neq__(self, nsphere):
return not self.__eq__(nsphere)
# ----------------------------------------------------------------------- #
# Size Setters/Getters #
# ----------------------------------------------------------------------- #
@property
def radius(self):
"""float: radius of n-sphere."""
return self.r
@property
def d(self):
"""float: diameter of n-sphere."""
return 2 * self.r
@property
def diameter(self):
"""float: diameter of n-sphere."""
return 2 * self.r
@property
def size(self):
"""float: size (diameter) of n-sphere."""
return 2 * self.r
@property
def position(self):
"""list: position of n-sphere."""
return self.center
# ----------------------------------------------------------------------- #
# Bounding N-Spheres #
# ----------------------------------------------------------------------- #
@property
def bound_max(self):
"""tuple: maximum bounding n-sphere"""
return tuple(list(self.center) + [self.r])
@property
def bound_min(self):
"""tuple: minimum interior n-sphere"""
return self.bound_max
# ----------------------------------------------------------------------- #
# Limits #
# ----------------------------------------------------------------------- #
@property
def limits(self):
"""list: list of (lower, upper) bounds for the bounding box"""
return [(x - self.r, x + self.r) for x in self.center]
@property
def sample_limits(self):
"""list: list of (lower, upper) bounds for the sampling region"""
return self.limits
# ----------------------------------------------------------------------- #
# Approximate #
# ----------------------------------------------------------------------- #
[docs] def approximate(self):
"""Approximate the n-sphere with itself
Other geometries can be approximated by a set of circles or spheres.
For the n-sphere, this approximation is exact.
Returns:
list: A list containing [(x, y, z, ..., r)]
"""
return [tuple(list(self.center) + [self.r])]
# ----------------------------------------------------------------------- #
# Within Test #
# ----------------------------------------------------------------------- #
[docs] def within(self, points):
"""Test if points are within n-sphere.
This function tests whether a point or set of points are within the
n-sphere. For the set of points, a list of booleans is returned to
indicate which points are within the n-sphere.
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)
sq_dist = np.sum(rel_pos * rel_pos, axis=-1)
mask = sq_dist <= self.r * self.r
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 n-sphere. Points at the center of the n-sphere 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)
cen_dist = np.sqrt(np.sum(rel_pos * rel_pos, axis=-1))
mask = cen_dist > 0
new_dist = 2 * self.r - cen_dist[mask]
scl = new_dist / cen_dist[mask]
new_rel_pos = np.zeros(pts.shape)
new_rel_pos[mask] = scl.reshape(-1, 1) * rel_pos[mask]
new_rel_pos[~mask] = 0
new_rel_pos[~mask, 0] = 2 * self.r
new_pts = new_rel_pos + np.array(self.center)
if single_pt:
return new_pts[0]
else:
return new_pts