# --------------------------------------------------------------------------- #
# #
# Import Modules #
# #
# --------------------------------------------------------------------------- #
from __future__ import division
import numpy as np
from matplotlib import patches
from matplotlib import pyplot as plt
from microstructpy.geometry.n_box import NBox
__author__ = 'Kenneth (Kip) Hart'
# --------------------------------------------------------------------------- #
# #
# Rectangle Class #
# #
# --------------------------------------------------------------------------- #
[docs]class Rectangle(NBox):
"""Rectangle
This class contains a generic, 2D rectangle. The position and dimensions
of the box can be specified using any of the parameters below.
Without parameters, this returns a unit square centered on the origin.
Args:
length (float): *(optional)* Length of the rectangle.
width (float): *(optional)* Width of the rectangle. *(optional)*
side_lengths (list): *(optional)* Side lengths. Defaults to (1, 1).
center (list): *(optional)* Center of rectangle. Defaults to (0, 0).
corner (list): *(optional)* Bottom-left corner.
bounds (list): *(optional)* Bounds of rectangle. Expected to be in the
format [(xmin, xmax), (ymin, ymax)].
limits : Alias for *bounds*.
angle (float): *(optional)* The rotation angle, in degrees.
angle_deg (float): *(optional)* The rotation angle, in degrees.
angle_rad (float): *(optional)* The rotation angle, in radians.
matrix (numpy.ndarray): *(optional)* The 2x2 rotation matrix.
"""
def __init__(self, **kwargs):
if 'length' in kwargs and 'width' in kwargs:
kwargs['side_lengths'] = [kwargs['length'], kwargs['width']]
if 'angle' in kwargs:
cp = np.cos(np.radians(kwargs['angle']))
sp = np.sin(np.radians(kwargs['angle']))
kwargs['matrix'] = np.array([[cp, -sp], [sp, cp]])
elif 'angle_deg' in kwargs:
cp = np.cos(np.radians(kwargs['angle_deg']))
sp = np.sin(np.radians(kwargs['angle_deg']))
kwargs['matrix'] = np.array([[cp, -sp], [sp, cp]])
elif 'angle_rad' in kwargs:
cp = np.cos(kwargs['angle_rad'])
sp = np.sin(kwargs['angle_rad'])
kwargs['matrix'] = np.array([[cp, -sp], [sp, cp]])
NBox.__init__(self, **kwargs)
try:
self.center
except AttributeError:
self.center = [0, 0]
try:
self.side_lengths
except AttributeError:
self.side_lengths = [1, 1]
# ----------------------------------------------------------------------- #
# Best Fit Function #
# ----------------------------------------------------------------------- #
[docs] def best_fit(self, points):
"""Find rectangle of best fit for points
Args:
points (list): List of points to fit.
Returns:
Rectangle: an instance of the class that best fits the points.
"""
# Unpack the input points
pts = np.array(points, dtype='float')
x, y = pts.T
# Find the most likely orientation for the rectangle
A = np.vstack([x, np.ones(len(x))]).T
m, _ = np.linalg.lstsq(A, y, rcond=None)[0]
theta = np.arctan(m)
# Rotate the points to an axis-aligned frame
s = np.sin(theta)
c = np.cos(theta)
rot = np.array([[c, -s], [s, c]])
pts_prin = pts.dot(rot)
# Translate points to center of bounding box
mins = pts_prin.min(axis=0)
maxs = pts_prin.max(axis=0)
bb_cen = 0.5 * (mins + maxs)
pts_bb = pts_prin - bb_cen
x_min, y_min = pts_bb.min(axis=0)
x_max, y_max = pts_bb.max(axis=0)
# Find nearest edge
x1_dist = np.abs(pts_bb[:, 0] - x_min)
x2_dist = np.abs(pts_bb[:, 0] - x_max)
y1_dist = np.abs(pts_bb[:, 1] - y_min)
y2_dist = np.abs(pts_bb[:, 1] - y_max)
dists = np.array([x1_dist, x2_dist, y1_dist, y2_dist]).T
min_dist = dists.min(axis=1)
mask_x1 = np.isclose(x1_dist, min_dist)
mask_x2 = np.isclose(x2_dist, min_dist)
mask_y1 = np.isclose(y1_dist, min_dist)
mask_y2 = np.isclose(y2_dist, min_dist)
# Group points by edge
x_x1 = pts_bb[mask_x1, 0].T
y_y1 = pts_bb[mask_y1, 1].T
x_x2 = pts_bb[mask_x2, 0].T
y_y2 = pts_bb[mask_y2, 1].T
# Get lines of best fit for each edge
x1_fit = np.mean(x_x1)
y1_fit = np.mean(y_y1)
x2_fit = np.mean(x_x2)
y2_fit = np.mean(y_y2)
x_cen = 0.5 * (x1_fit + x2_fit)
y_cen = 0.5 * (y1_fit + y2_fit)
fit_cen_bb = np.array([x_cen, y_cen])
x_len = x2_fit - x1_fit
y_len = y2_fit - y1_fit
# Translate center to rotated frame
fit_cen_prin = fit_cen_bb + bb_cen
fit_cen = rot.dot(fit_cen_prin.reshape(-1, 1)).flatten()
# Disambiguate the orientation and axes
x_ax_seed = np.array(self.matrix)[:, 0]
x_dot, y_dot = rot.T.dot(x_ax_seed)
if np.abs(x_dot) > np.abs(y_dot):
if x_dot > 0:
x_ax_fit = rot[:, 0]
else:
x_ax_fit = - rot[:, 0]
length = x_len
width = y_len
else:
if y_dot > 0:
x_ax_fit = rot[:, 1]
else:
x_ax_fit = - rot[:, 1]
length = y_len
width = x_len
ang_diff = np.arcsin(np.cross(x_ax_seed, x_ax_fit))
ang_rad = self.angle_rad + ang_diff
rect_fit = Rectangle(center=fit_cen, length=length, width=width,
angle_rad=ang_rad)
return rect_fit
# ----------------------------------------------------------------------- #
# Representation Function #
# ----------------------------------------------------------------------- #
def __repr__(self):
repr_str = 'Rectangle('
repr_str += 'center=' + repr(tuple(self.center)) + ', '
repr_str += 'side_lengths=' + repr(tuple(self.side_lengths)) + ', '
repr_str += 'angle=' + repr(self.angle) + ')'
return repr_str
# ----------------------------------------------------------------------- #
# Number of Dimensions #
# ----------------------------------------------------------------------- #
@property
def n_dim(self):
"""int: Number of dimensions, 2"""
return 2
# ----------------------------------------------------------------------- #
# Area #
# ----------------------------------------------------------------------- #
@property
def area(self):
"""float: Area of rectangle"""
return self.n_vol
[docs] @classmethod
def area_expectation(cls, **kwargs):
r"""Expected area of rectangle
This method computes the expected area of a rectangle. There are two
main ways to define the size of a rectangle: by the length and width
and by the bounds. If the input rectangle is defined by length and
width, the expected area is:
.. math::
\mathbb{E}[A] = \mathbb{E}[L W] = \mu_L \mu_W
For the case where it is defined by upper and lower bounds:
.. math::
\mathbb{E}[A] = \mathbb{E}[(X_U - X_L) (Y_U - Y_L)]
.. math::
\mathbb{E}[A] =
\mu_{X_U}\mu_{Y_U} - \mu_{X_U} \mu_{Y_L} -
\mu_{X_L}\mu_{Y_U} + \mu_{X_L}\mu_{Y_L}
Example:
>>> import scipy.stats
>>> import microstructpy as msp
>>> L = scipy.stats.uniform(loc=1, scale=2)
>>> W = scipy.stats.norm(loc=3.2, scale=5.1)
>>> L.mean() * W.mean()
6.4
>>> msp.geometry.Rectangle.area_expectation(length=L, width=W)
6.4
Args:
**kwargs: Keyword arguments, same as :class:`.Rectangle` but the
inputs can be from the :mod:`scipy.stats` module.
Returns:
float: Expected/average area of rectangle.
"""
if 'length' in kwargs or 'width' in kwargs:
len_dist = kwargs.get('length', 1)
width_dist = kwargs.get('width', 1)
return _prod_exp(*[len_dist, width_dist])
if 'side_lengths' in kwargs:
return _prod_exp(*kwargs['side_lengths'])
if 'bounds' in kwargs:
x_bnds, y_bnds = kwargs['bounds']
x_lb, x_ub = x_bnds
y_lb, y_ub = y_bnds
p1 = _prod_exp(*[x_ub, y_ub])
p2 = _prod_exp(*[x_ub, y_lb])
p3 = _prod_exp(*[x_lb, y_ub])
p4 = _prod_exp(*[x_lb, y_lb])
return p1 - p2 - p3 + p4
raise ValueError('Cannot find the expected area of rectangle')
# ----------------------------------------------------------------------- #
# Properties #
# ----------------------------------------------------------------------- #
@property
def length(self):
"""float: Length of rectangle, side length along 1st axis"""
return self.side_lengths[0]
@property
def width(self):
"""float: Width of rectangle, side length along 2nd axis"""
return self.side_lengths[1]
@property
def angle(self):
"""float: Rotation angle of rectangle - degrees"""
return self.angle_deg
@property
def angle_deg(self):
"""float: Rotation angle of rectangle - degrees"""
return 180 * self.angle_rad / np.pi
@property
def angle_rad(self):
"""float: Rotation angle of rectangle - radians"""
return np.arctan2(self.matrix[1][0], self.matrix[0][0])
# ----------------------------------------------------------------------- #
# Circle Approximation #
# ----------------------------------------------------------------------- #
[docs] def approximate(self, x1=None):
"""Approximate rectangle with a set of circles.
This method approximates a rectangle with a set of circles.
These circles are spaced uniformly along the long axis of the
rectangle with distance ``x1`` between them.
Example
-------
For a rectangle with length=2.5, width=1, and x1=0.3,
the approximation would look like :numref:`f_api_rectangle_approx`.
.. _f_api_rectangle_approx:
.. figure:: ../../auto_examples/geometry/images/sphx_glr_plot_rectangle_001.png
Circular approximation of rectangle.
Args:
x1 (float or None): *(optional)* Spacing between the circles.
If not specified, the spacing is 0.25x the length of the
shortest side.
Returns:
numpy.ndarray: An Nx3 array, where each row is a circle and the
columns are x, y, and r.
""" # NOQA: E501
if x1 is None:
x1 = 0.25 * min(self.side_lengths)
if self.side_lengths[0] >= self.side_lengths[1]:
length = self.side_lengths[0]
width = self.side_lengths[1]
inds = [0, 1]
else:
length = self.side_lengths[1]
width = self.side_lengths[0]
inds = [1, 0]
# Centerline circles
xc = 0
half_len = 0.5 * length
r = 0.5 * width
circs = []
while xc < half_len:
circ = [xc, 0, r]
circs.append(circ)
if np.isclose(xc + r, half_len):
xc = half_len
elif xc + x1 > half_len - r:
xc = half_len - r
else:
xc = xc + x1
# Corner circle
while r > 0:
x_init = circs[-1][0]
y_init = circs[-1][1]
x = x_init + x1
y = y_init + x1
dx = half_len - x
dy = 0.5 * width - y
r = min(dx, dy)
if r <= 0:
break
else:
circs.append([x, y, r])
# Reflect circles
circs = np.array(circs)
for dim in range(self.n_dim):
mask = circs[:, dim] > 0
new_circs = np.copy(circs[mask])
new_circs[:, dim] *= -1
circs = np.concatenate((circs, new_circs))
circs[:, inds] = circs[:, [0, 1]]
# Rotate and translate circles
pts = circs[:, :-1]
circs[:, :-1] = pts.dot(np.array(self.matrix).T)
circs[:, :-1] += self.center
return circs
# ----------------------------------------------------------------------- #
# Plot Function #
# ----------------------------------------------------------------------- #
[docs] def plot(self, **kwargs):
"""Plot the rectangle.
This function adds a :class:`matplotlib.patches.Rectangle` patch to the
current axes. The keyword arguments are passed through to the patch.
Args:
**kwargs (dict): Keyword arguments for the patch.
""" # NOQA: E501
pt = self.corner
w = self.length
h = self.width
ang = self.angle
c = patches.Rectangle(pt, w, h, angle=ang, **kwargs)
plt.gca().add_patch(c)
def _prod_exp(*args):
prod = 1
for arg in args:
try:
arg_mu = arg.moment(1)
except AttributeError:
arg_mu = arg
prod *= arg_mu
return prod
# --------------------------------------------------------------------------- #
# #
# Square Class #
# #
# --------------------------------------------------------------------------- #
[docs]class Square(Rectangle):
"""A square.
This class contains a generic, 2D square. It is derived from the
:class:`microstructpy.geometry.Rectangle` class and contains the
``side_length`` property, rather than multiple side lengths.
Args:
side_length (float): *(optional)* Side length. Defaults to 1.
center (list): *(optional)* Center of rectangle. Defaults to (0, 0).
corner (list): *(optional)* Bottom-left corner.
"""
def __init__(self, **kwargs):
if 'side_length' in kwargs:
kwargs['side_lengths'] = 2 * [kwargs['side_length']]
Rectangle.__init__(self, **kwargs)
# ----------------------------------------------------------------------- #
# Side Length Property #
# ----------------------------------------------------------------------- #
@property
def side_length(self):
"""float: length of the side of the square."""
return self.side_lengths[0]
# ----------------------------------------------------------------------- #
# Area #
# ----------------------------------------------------------------------- #
[docs] @classmethod
def area_expectation(cls, **kwargs):
r"""Expected area of square
This method computes the expected area of a square with distributed
side length.
The expectation is:
.. math::
\mathbb{E}[A] = \mathbb{E}[S^2] = \mu_S^2 + \sigma_S^2
Example:
>>> import scipy.stats
>>> import microstructpy as msp
>>> S = scipy.stats.expon(scale=2)
>>> S.mean()^2 + S.var()
8.0
>>> msp.geometry.Square.area_expectation(side_length=S)
8.0
Args:
**kwargs: Keyword arguments, same as :class:`.Square` but the
inputs can be from the :mod:`scipy.stats` module.
Returns:
float: Expected/average area of the square.
"""
if 'side_length' in kwargs:
len_dist = kwargs['side_length']
try:
area_exp = len_dist.moment(2)
except AttributeError:
area_exp = len_dist * len_dist
return area_exp
Rectangle.area_expectation(**kwargs)
# ----------------------------------------------------------------------- #
# Circle Approximation #
# ----------------------------------------------------------------------- #
[docs] def approximate(self, x1=None):
"""Approximate square with a set of circles
This method approximates a square with a set of circles.
These circles are spaced uniformly along the edges of the square
with distance ``x1`` between them.
Example
-------
For a square with side_length=1, and x1=0.2,
the approximation would look like :numref:`f_api_square_approx`.
.. _f_api_square_approx:
.. figure:: ../../auto_examples/geometry/images/sphx_glr_plot_rectangle_002.png
Circular approximation of square.
Args:
x1 (float or None): *(optional)* Spacing between the circles.
If not specified, the spacing is 0.25x the side length.
Returns:
numpy.ndarray: An Nx3 array, where each row is a circle and the
columns are x, y, and r.
""" # NOQA: E501
return Rectangle.approximate(self, x1)
