# Source code for sknano.core.geometric_regions._3D_regions

# -*- coding: utf-8 -*-
"""
=======================================================================
3D geometric regions (:mod:sknano.core.geometric_regions._3D_regions)
=======================================================================

.. currentmodule:: sknano.core.geometric_regions._3D_regions

"""
from __future__ import absolute_import, division, print_function, \
unicode_literals
__docformat__ = 'restructuredtext en'

from abc import ABCMeta, abstractmethod

import numbers
import numpy as np

from sknano.core.math import Point, Vector, vector as vec
from ._base import GeometricRegion, GeometricTransformsMixin

__all__ = ['Geometric3DRegion', 'Parallelepiped', 'Cuboid', 'Cube',
'Ellipsoid', 'Sphere', 'Cylinder', 'Cone']

[docs]class Geometric3DRegion(GeometricRegion, GeometricTransformsMixin,
metaclass=ABCMeta):
"""Abstract base class for representing 3D geometric regions."""

@property
@abstractmethod
def volume(self):
"""Volume of 3D geometric region."""
raise NotImplementedError

@property
def measure(self):
"""Alias for :attr:~Geometric3DRegion.volume, which is the measure \
of a 3D geometric region."""
return self.volume

[docs]class Parallelepiped(Geometric3DRegion):
"""Geometric3DRegion for a parallelepiped.

Represents a parallelepiped with origin :math:(o_x, o_y, o_z) and
direction vectors :math:\\mathbf{u}=(u_x, u_y, u_z),
:math:\\mathbf{v}=(v_x, v_y, v_z), and
:math:\\mathbf{w}=(w_x, w_y, w_z).

Parameters
----------
o : array_like, optional
Parallelepiped origin. If None, it defaults to o=[0, 0, 0].
u, v, w : array_like, optional
Parallelepiped direction vectors stemming from origin o.
If None, then the default values are u=[1, 0, 0],
v=[1, 1, 0], and w=[0, 1, 1].

Notes
-----
:class:Parallelepiped represents the bounded region
:math:\\left\\{o+\\lambda_1\\mathbf{u}+\\lambda_2\\mathbf{v}+
\\lambda_3\\mathbf{w}\\in R^3|0\\le\\lambda_i\\le 1\\right\\},
where :math:\\mathbf{u}, :math:\\mathbf{v}, and :math:\\mathbf{w}
have to be linearly independent.

Calling Parallelepiped with no parameters is equivalent to
Parallelepiped\ (o=[0, 0, 0], u=[1, 0, 0], v=[1, 1, 0], w=[0, 1, 1]).

"""
def __init__(self, o=None, u=None, v=None, w=None):

self._o = Point()
self._u = Vector()
self._v = Vector()
self._w = Vector()

super().__init__()

if o is None:
o = [0, 0, 0]
self.o = o

if u is None:
u = [1, 0, 0]
self.u = u

if v is None:
v = [1, 1, 0]
self.v = v

if w is None:
w = [0, 1, 1]
self.w = w

self.points.append(self.o)
self.vectors.extend([self.u, self.v, self.w])

self.fmtstr = "o={o!r}, u={u!r}, v={v!r}, w={w!r}"

@property
def o(self):
"""3D point coordinates :math:(o_x, o_y, o_z) of origin.

Returns
-------
:class:Point
3D :class:Point coordinates :math:(o_x, o_y, o_z) of origin.
"""
return self._o

@o.setter
def o(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self.vectors.translate(Vector(p0=self.o, p=Point(value)))
self._o[:] = Point(value)

@property
def u(self):
"""3D direction vector :math:\\mathbf{u}=(u_x, u_y, u_z), with \
origin :attr:~Parallelepiped.o"""
return self._u

@u.setter
def u(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._u[:] = Vector(value, p0=self.o)

@property
def v(self):
"""3D direction vector :math:\\mathbf{v}=(v_x, v_y, v_z), with \
origin :attr:~Parallelepiped.o"""
return self._v

@v.setter
def v(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._v[:] = Vector(value, p0=self.o)

@property
def w(self):
"""3D direction vector :math:\\mathbf{w}=(w_x, w_y, w_z), with \
origin :attr:~Parallelepiped.o"""
return self._w

@w.setter
def w(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._w[:] = Vector(value, p0=self.o)

@property
def volume(self):
"""Parallelepiped volume, :math:V.

Computed as:

.. math::

V = |\\mathbf{u}\\cdot\\mathbf{v}\\times\\mathbf{w}|

"""
return np.abs(vec.scalar_triple_product(self.u, self.v, self.w))

@property
def centroid(self):
"""Parallelepiped centroid, :math:(c_x, c_y, c_z).

Computed as the 3D point :math:(c_x, c_y, c_z) with coordinates:

.. math::

c_x = o_x + \\frac{u_x + v_x + w_x}{2}

c_y = o_y + \\frac{u_y + v_y + w_y}{2}

c_z = o_z + \\frac{u_z + v_z + w_z}{2}

where :math:(o_x, o_y, o_z), :math:(u_x, u_y, u_z),
:math:(v_x, v_y, v_z), and :math:(w_x, w_y, w_z)
are the :math:(x, y, z) coordinates of the origin :math:o
and :math:(x, y, z) components of the direction vectors
:math:\\mathbf{u}, :math:\\mathbf{v}, and :math:\\mathbf{w},
respectively.

Returns
-------
:class:Point
3D :class:Point of centroid.

"""
ox, oy, oz = self.o
ux, uy, uz = self.u
vx, vy, vz = self.v
wx, wy, wz = self.w

cx = ox + (ux + vx + wx) / 2
cy = oy + (uy + vy + wy) / 2
cz = oz + (uz + vz + wz) / 2

return Point([cx, cy, cz])

[docs]    def contains(self, point):
"""Test region membership of point in :class:Parallelepiped.

Parameters
----------
point : array_like

Returns
-------
:class:~python:bool
True if point is within :class:Parallelepiped,
False otherwise

Notes
-----
A point :math:(p_x, p_y, p_z) is within the bounded region
of a parallelepiped with origin :math:(o_x, o_y, o_z) and
direction vectors :math:\\mathbf{u}=(u_x, u_y, u_z),
:math:\\mathbf{v}=(v_x, v_y, v_z), and
:math:\\mathbf{w}=(w_x, w_y, w_z) if the following is true:

.. math::

0\\le\\frac{
v_z (w_x (p_y - o_y) + w_y (o_x - p_x)) +
w_z (v_x (o_y - p_y) + v_y (p_x - o_x)) +
o_z (v_y w_x - v_x w_y) + p_z (v_x w_y - v_y w_x)}{
u_z (v_x w_y - v_y w_x) +
u_y (v_z w_x - v_x w_z) +
u_x (v_y w_z - v_z w_y)}\\le 1 \\land

0\\le\\frac{
u_z (w_x (p_y - o_y) + w_y (o_x - p_x)) +
w_z (u_x (o_y - p_y) + u_y (p_x - o_x)) +
o_z (u_y w_x - u_x w_y) +
p_z (u_x w_y - u_y w_x)}{
u_z (v_y w_x - v_x w_y) +
u_y (v_x w_z - v_z w_x) +
u_x (v_z w_y - v_y w_z)}\\le 1 \\land

0\\le\\frac{
u_z (v_x (p_y - o_y) + v_y (o_x - p_x)) +
v_z (u_x (o_y - p_y) + u_y (p_x - o_x)) +
o_z (u_y v_x - u_x v_y) +
p_z (u_x v_y - u_y v_x)}{
u_z (v_x w_y - v_y w_x) +
u_y (v_z w_x - v_x w_z) +
u_x (v_y w_z - v_z w_y)}\\le 1

"""
px, py, pz = Point(point)

ox, oy, oz = self.o
ux, uy, uz = self.u
vx, vy, vz = self.v
wx, wy, wz = self.w

q1 = (vz * (wx * (py - oy) + wy * (ox - px)) +
wz * (vx * (oy - py) + vy * (px - ox)) +
oz * (vy * wx - vx * wy) +
pz * (vx * wy - vy * wx)) / \
(uz * (vx * wy - vy * wx) +
uy * (vz * wx - vx * wz) +
ux * (vy * wz - vz * wy))

q2 = (uz * (wx * (py - oy) + wy * (ox - px)) +
wz * (ux * (oy - py) + uy * (px - ox)) +
oz * (uy * wx - ux * wy) +
pz * (ux * wy - uy * wx)) / \
(uz * (vy * wx - vx * wy) +
uy * (vx * wz - vz * wx) +
ux * (vz * wy - vy * wz))

q3 = (uz * (vx * (py - oy) + vy * (ox - px)) +
vz * (ux * (oy - py) + uy * (px - ox)) +
oz * (uy * vx - ux * vy) +
pz * (ux * vy - uy * vx)) / \
(uz * (vx * wy - vy * wx) +
uy * (vz * wx - vx * wz) +
ux * (vy * wz - vz * wy))

return q1 >= 0 and q1 <= 1 and q2 >= 0 and q2 <= 1 and \
q3 >= 0 and q3 <= 1

[docs]    def todict(self):
"""Returns a :class:~python:dict of the :class:Parallelepiped \
constructor parameters."""
return dict(o=self.o, u=self.u, v=self.v, w=self.w)

[docs]class Cuboid(Geometric3DRegion):
"""Geometric3DRegion for a cuboid.

Represents an axis-aligned cuboid with lower corner
:math:p_{\\mathrm{min}}=
(x_{\\mathrm{min}},y_{\\mathrm{min}},z_{\\mathrm{min}}) and
upper corner
:math:p_{\\mathrm{max}}=
(x_{\\mathrm{max}},y_{\\mathrm{max}},z_{\\mathrm{max}}).

Parameters
----------
pmin, pmax : array_like, optional
The minimum and maximum 3D point coordinates of the axis-aligned
cuboid from pmin=[xmin, ymin, zmin] to pmax=[xmax, ymax, zmax].
xmin, ymin, zmin : float, optional
The minimum :math:(x, y, z) coordinates of the axis-aligned cuboid.
xmax, ymax, zmax : float, optional
The maximum :math:(x, y, z) coordinates of the axis-aligned cuboid.

Notes
-----
:class:Cuboid represents the region
:math:\\left\\{\\{x, y, z\\}|
x_{\\mathrm{min}}\\le x\\le x_{\\mathrm{max}}\\land
y_{\\mathrm{min}}\\le y\\le y_{\\mathrm{max}}\\land
z_{\\mathrm{min}}\\le z\\le z_{\\mathrm{max}}\\right\\}

Calling Cuboid with no parameters is equivalent to
Cuboid\ (pmin=[0, 0, 0], pmax=[1, 1, 1]).

"""
def __init__(self, pmin=None, pmax=None, xmin=0, ymin=0, zmin=0,
xmax=1, ymax=1, zmax=1):
super().__init__()
self._pmin = Point()
self._pmax = Point()

if pmin is None:
pmin = [xmin, ymin, zmin]
self.pmin = pmin

if pmax is None:
pmax = [xmax, ymax, zmax]
self.pmax = pmax

# assert self.pmin <= self.pmax
assert np.all(np.less_equal(self.pmin, self.pmax))

self.points.extend([self.pmin, self.pmax])
self.fmtstr = "pmin={pmin!r}, pmax={pmax!r}"

@property
def pmin(self):
"""3D :class:Point at \
(:attr:~Cuboid.xmin, :attr:~Cuboid.ymin, :attr:~Cuboid.zmin)
"""
return self._pmin

@pmin.setter
def pmin(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._pmin[:] = Point(value)

@property
def pmax(self):
"""3D :class:Point at \
(:attr:~Cuboid.xmax, :attr:~Cuboid.ymax, :attr:~Cuboid.zmax)
"""
return self._pmax

@pmax.setter
def pmax(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._pmax[:] = Point(value)

@property
def xmin(self):
""":math:x_{\\mathrm{min}} coordinate."""
return self.pmin.x

@xmin.setter
def xmin(self, value):
if not isinstance(value, numbers.Number):
raise TypeError('Expected a number')
self.pmin.x = float(value)

@property
def ymin(self):
""":math:y_{\\mathrm{min}} coordinate."""
return self.pmin.y

@ymin.setter
def ymin(self, value):
if not isinstance(value, numbers.Number):
raise TypeError('Expected a number')
self.pmin.y = float(value)

@property
def zmin(self):
""":math:z_{\\mathrm{min}} coordinate."""
return self.pmin.z

@zmin.setter
def zmin(self, value):
if not isinstance(value, numbers.Number):
raise TypeError('Expected a number')
self.pmin.z = float(value)

@property
def xmax(self):
""":math:x_{\\mathrm{max}} coordinate."""
return self.pmax.x

@xmax.setter
def xmax(self, value):
if not isinstance(value, numbers.Number):
raise TypeError('Expected a number')
self.pmax.x = float(value)

@property
def ymax(self):
""":math:y_{\\mathrm{max}} coordinate."""
return self.pmax.y

@ymax.setter
def ymax(self, value):
if not isinstance(value, numbers.Number):
raise TypeError('Expected a number')
self.pmax.y = float(value)

@property
def zmax(self):
""":math:z_{\\mathrm{max}} coordinate."""
return self.pmax.z

@zmax.setter
def zmax(self, value):
if not isinstance(value, numbers.Number):
raise TypeError('Expected a number')
self.pmax.z = float(value)

@property
def a(self):
"""Distance between :math:x_{\\mathrm{max}}-x_{\\mathrm{min}}."""
return self.xmax - self.xmin

@property
def b(self):
"""Distance between :math:y_{\\mathrm{max}}-y_{\\mathrm{min}}."""
return self.ymax - self.ymin

@property
def c(self):
"""Distance between :math:z_{\\mathrm{max}}-z_{\\mathrm{min}}."""
return self.zmax - self.zmin

@property
def volume(self):
""":class:Cuboid volume, :math:V=abc."""
return self.a * self.b * self.c

@property
def centroid(self):
""":class:Cuboid centroid, :math:(c_x, c_y, c_z).

Computed as the 3D :class:Point :math:(c_x, c_y, c_z) with
coordinates:

.. math::

c_x = \\frac{x_{\\mathrm{min}} + x_{\\mathrm{max}}}{2}

c_y = \\frac{y_{\\mathrm{min}} + y_{\\mathrm{max}}}{2}

c_z = \\frac{z_{\\mathrm{min}} + z_{\\mathrm{max}}}{2}

Returns
-------
:class:Point
3D :class:Point of centroid.

"""
cx = (self.xmin + self.xmax) / 2
cy = (self.ymin + self.ymax) / 2
cz = (self.zmin + self.zmax) / 2
return Point([cx, cy, cz])

[docs]    def contains(self, point):
"""Test region membership of point in :class:Cuboid.

Parameters
----------
point : array_like

Returns
-------
:class:~python:bool
True if point is within :class:Cuboid,
False, otherwise.

Notes
-----
A point :math:(p_x, p_y, p_z) is within the bounded region of a
cuboid with lower corner at
:math:p_{\\mathrm{min}}=
(x_{\\mathrm{min}}, y_{\\mathrm{min}}, z_{\\mathrm{min}}) and
upper corner at
:math:p_{\\mathrm{max}}=
(x_{\\mathrm{max}}, y_{\\mathrm{max}}, z_{\\mathrm{max}}) if the
following is true:

.. math::

x_{\mathrm{min}}\\le x\\le x_{\\mathrm{max}}\\land

y_{\mathrm{min}}\\le y\\le y_{\\mathrm{max}}\\land

z_{\mathrm{min}}\\le z\\le z_{\\mathrm{max}}

"""
px, py, pz = Point(point)

return (px >= self.xmin) and (px <= self.xmax) and \
(py >= self.ymin) and (py <= self.ymax) and \
(pz >= self.zmin) and (pz <= self.zmax)

[docs]    def rotate(self, **kwargs):
super().rotate(**kwargs)
pmin = self.pmin.copy()
pmax = self.pmax.copy()
self.pmin = np.minimum(pmin, pmax)
self.pmax = np.maximum(pmin, pmax)
assert np.all(np.less_equal(self.pmin, self.pmax))

[docs]    def todict(self):
"""Returns a :class:~python:dict of the :class:Cuboid \
constructor parameters."""
return dict(pmin=self.pmin, pmax=self.pmax)

[docs]class Cube(Geometric3DRegion):
"""Geometric3DRegion for a cube.

:class:Cube represents the region
:math:\\left\\{\\{c_i\\pm\\frac{a}{2}\\}|a>0\\forall
i\\in\\{x,y,z\\}\\right\\}.

Parameters
----------
center : array_like, optional
The :math:(x,y,z) coordinate of the axis-aligned cube
center :math:(c_x, c_y, c_z).
a : float, optional
Side length :math:a of axis-aligned cube.

Notes
-----
Calling :class:Cube with no parameters is equivalent to
:class:Cube\ (center=[0, 0, 0], a=1).

"""
def __init__(self, center=None, a=1):
super().__init__()

self._center = Point()

if center is None:
center = [0, 0, 0]
self.center = center
self.a = a
self.points.append(self.center)

self.fmtstr = "center={center!r}, a={a:.3f}"

@property
def center(self):
"""Center point :math:(c_x, c_y, c_z) of axis-aligned cube."""
return self._center

@center.setter
def center(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._center[:] = Point(value)

@property
def a(self):
"""Side length :math:a of the axis-aligned cube."""
return self._a

@a.setter
def a(self, value):
self._a = float(value)

@property
def centroid(self):
"""Alias for :attr:~Cube.center."""
return self.center

@property
def volume(self):
""":class:Cube volume, :math:V=a^3."""
return self.a ** 3

[docs]    def contains(self, point):
"""Test region membership of point in :class:Cube.

Parameters
----------
point : array_like

Returns
-------
:class:~python:bool
True if point is within :class:Cube,
False otherwise

Notes
-----
A point :math:(p_x, p_y, p_z) is within the bounded region of
a cube with center :math:(c_x, c_y, c_z) and side length
:math:a if the following is true:

.. math::

c_i-\\frac{a}{2}\\le p_i\\le c_i+\\frac{a}{2}\\forall
i\\in \\{x, y, z\\}

"""
px, py, pz = Point(point)
cx, cy, cz = self.center
a = self.a
xmin = cx - a / 2
xmax = cx + a / 2
ymin = cy - a / 2
ymax = cy + a / 2
zmin = cz - a / 2
zmax = cz + a / 2
return (px >= xmin) and (px <= xmax) and \
(py >= ymin) and (py <= ymax) and \
(pz >= zmin) and (pz <= zmax)

[docs]    def todict(self):
"""Returns a :class:~python:dict of the :class:Cube \
constructor parameters."""
return dict(center=self.center, a=self.a)

[docs]class Ellipsoid(Geometric3DRegion):
"""Geometric3DRegion for an ellipsoid.

Represents an axis-aligned ellipsoid centered at the point
:math:(c_x, c_y, c_z) with semi-axes lengths
:math:r_x, r_y, r_z.

Parameters
----------
center : array_like or :class:~sknano.core.Point
3-tuple of floats or an instance of the
:class:~sknano.core.Point class specifying the center point
coordinates :math:(x,y,z) of the axis-aligned
:class:Ellipsoid with semi-axes lengths rx, ry, rz.
rx, ry, rz : float, optional
Semi-axes lengths :math:r_x, r_y, r_z of axis-aligned
:class:Ellipsoid centered at center.

Notes
-----
:class:Ellipsoid represents the axis-aligned ellipsoid:

.. math::

\\left\\{\\{x, y, z\\}\\in R^3|
\\left(\\frac{x-c_x}{r_x}\\right)^2 +
\\left(\\frac{y-c_y}{r_y}\\right)^2 +
\\left(\\frac{z-c_z}{r_z}\\right)^2\\le 1\\right\\}

Calling :class:Ellipsoid with no parameters is equivalent to
:class:Ellipsoid\ (center=[0, 0, 0], rx=1, ry=1, rz=1).

"""
def __init__(self, center=None, rx=1, ry=1, rz=1):
super().__init__()

self._center = Point()

if center is None:
center = [0, 0, 0]
self.center = center

self.rx = rx
self.ry = ry
self.rz = rz

self.points.append(self.center)

self.fmtstr = \
"center={center!r}, rx={rx:.3f}, ry={ry:.3f}, rz={rz:.3f}"

@property
def center(self):
""":class:Ellipsoid center point :math:(c_x, c_y, c_z)."""
return self._center

@center.setter
def center(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._center[:] = Point(value)

@property
def rx(self):
"""Length of semi-axis :math:r_x."""
return self._rx

@rx.setter
def rx(self, value):
self._rx = float(value)

@property
def ry(self):
"""Length of semi-axis :math:r_y."""
return self._ry

@ry.setter
def ry(self, value):
self._ry = float(value)

@property
def rz(self):
"""Length of semi-axis :math:r_z."""
return self._rz

@rz.setter
def rz(self, value):
self._rz = float(value)

@property
def centroid(self):
"""Alias for :attr:~Ellipsoid.center."""
return self.center

@property
def volume(self):
""":class:Ellipsoid volume, \
:math:V=\\frac{4}{3}\\pi r_x r_y r_z."""
return 4 / 3 * np.pi * self.rx * self.ry * self.rz

[docs]    def contains(self, point):
"""Test region membership of point in :class:Ellipsoid.

Parameters
----------
point : array_like

Returns
-------
:class:~python:bool
True if point is within :class:Ellipsoid,
False otherwise

Notes
-----
A point :math:(p_x, p_y, p_z) is within the bounded region of
an ellipsoid with center :math:(c_x, c_y, c_z) and semi-axes lengths
:math:r_x, r_y, r_z if the following is true:

.. math::

\\left(\\frac{p_x - c_x}{r_x}\\right)^2 +
\\left(\\frac{p_y - c_y}{r_y}\\right)^2 +
\\left(\\frac{p_z - c_z}{r_z}\\right)^2\\le 1

"""
px, py, pz = Point(point)
cx, cy, cz = self.center
rx, ry, rz = self.rx, self.ry, self.rz

q1 = (px - cx) ** 2 / rx ** 2
q2 = (py - cy) ** 2 / ry ** 2
q3 = (pz - cz) ** 2 / rz ** 2

return q1 + q2 + q3 <= 1.0

[docs]    def todict(self):
"""Returns a :class:~python:dict of the :class:Ellipsoid \
constructor parameters."""
return dict(center=self.center, rx=self.rx, ry=self.ry, rz=self.rz)

[docs]class Sphere(Geometric3DRegion):
"""Geometric3DRegion for a sphere.

Parameters
----------
center : array_like, optional
Either a 3-tuple of floats or an instance of the
:class:~sknano.core.Point class specifying the :math:(x,y,z)
coordinates of the :class:Sphere center.
r : float, optional
Sphere radius :math:r
"""
def __init__(self, center=None, r=1):
super().__init__()

self._center = Point()

if center is None:
center = [0, 0, 0]
self.center = center

self.r = r

self.points.append(self.center)

self.fmtstr = "center={center!r}, r={r:.3f}"

@property
def center(self):
""":class:Sphere center point :math:(h, k, l)."""
return self._center

@center.setter
def center(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._center[:] = Point(value)

@property
def centroid(self):
"""Alias for :attr:~Sphere.center."""
return self.center

@property
def r(self):
""":class:Sphere radius, :math:r."""
return self._r

@r.setter
def r(self, value):
self._r = float(value)

@property
def volume(self):
""":class:Sphere volume, :math:V=\\frac{4}{3}\\pi r^3."""
return 4 / 3 * np.pi * self.r ** 3

[docs]    def contains(self, point):
"""Test region membership of point in :class:Sphere.

Parameters
----------
point : array_like

Returns
-------
:class:~python:bool
True if point is within :class:Sphere,
False otherwise.

Notes
-----
A point :math:(p_x, p_y, p_z) is within the bounded region
of a sphere with center :math:(h, k, l) and radius :math:r
if the following is true:

.. math::

(p_x - h)^2 + (p_y - k)^2 + (p_z - l)^2 \\le r^2

"""
x, y, z = Point(point)
h, k, l = self.center
r = self.r

return (x - h) ** 2 + (y - k) ** 2 + (z - l) ** 2 <= r ** 2

[docs]    def todict(self):
"""Returns a :class:~python:dict of the :class:Sphere \
constructor parameters."""
return dict(center=self.center, r=self.r)

[docs]class Cylinder(Geometric3DRegion):
"""Geometric3DRegion for a cylinder.

Represents a cylinder of radius :math:r around the line from
:math:(x_1, y_1, z_1) to :math:(x_2, y_2, z_2).

Parameters
----------
p1, p2 : array_like, optional
3-tuples or :class:~sknano.core.Point class instances
specifying the Cylinder axis from point :math:p_1=(x_1,y_1,z_1)
to :math:p_2=(x_2,y_2,z_2).
r : float, optional
Cylinder radius :math:r

Notes
-----
:class:Cylinder represents a cylinder region
:math:\\left\\{p_1+\\rho\\cos(\\theta)\\mathbf{v}_1 +
\\rho\\sin(\\theta)\\mathbf{v}_2 + \\mathbf{v}_3 z|
0\\le\\theta\\le 2\\pi\\land
0\\le\\rho\\le 1\\land
0\\le z\\le 1\\right\\}
where :math:\\mathbf{v}_3=p_2 - p_1 and the vectors
:math:\\{\\mathbf{v}_1,\\mathbf{v}_2,\\mathbf{v}_3\\} are orthogonal
with :math:|\\mathbf{v}_1|=|\\mathbf{v}_2|=1, and
:math:p_1=(x_1, y_1, z_1) and :math:p_2=(x_2,y_2,z_2).

Calling :class:Cylinder with no parameters is equivalent to
:class:Cylinder\ (p1=[0, 0, -1], p2=[0, 0, 1], r=1).

"""
def __init__(self, p1=None, p2=None, r=1):
super().__init__()

self._p1 = Point()
self._p2 = Point()

if p1 is None:
p1 = [0, 0, -1]
if p2 is None:
p2 = [0, 0, 1]

self.p1 = p1
self.p2 = p2
self.points.extend([self.p1, self.p2])
self.r = r
self.fmtstr = "p1={p1!r}, p2={p2!r}, r={r:.3f}"

@property
def r(self):
""":class:Cylinder radius :math:r."""
return self._r

@r.setter
def r(self, value):
self._r = float(value)

@property
def p1(self):
""":class:Cylinder axis point :math:p_1=(x_1, y_1, z_1)."""
return self._p1

@p1.setter
def p1(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._p1[:] = Point(value)

@property
def p2(self):
""":class:Cylinder axis point :math:p_2=(x_2, y_2, z_2)."""
return self._p2

@p2.setter
def p2(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._p2[:] = Point(value)

@property
def axis(self):
""":class:Cylinder axis :class:Vector \
:math:\\boldsymbol{\\ell}=p_2 - p_1.

Returns
-------
:class:Vector
3D :class:Vector along :class:Cylinder axis from
:class:Point :math:p_1=(x_1, y_1, z_1) to
:class:Point :math:p_2=(x_2, y_2, z_2).

"""
return Vector(p0=self.p1, p=self.p2)

@property
def centroid(self):
""":class:Cylinder centroid, :math:(c_x, c_y, c_z).

Computed as:

.. math::

c_x = \\frac{x_1 + x_2}{2}

c_y = \\frac{y_1 + y_2}{2}

c_z = \\frac{z_1 + z_2}{2}

"""
return (self.p1 + self.p2) / 2

@property
def volume(self):
""":class:Cylinder volume, :math:V=\\pi r^2 \\ell."""
return np.pi * self.r ** 2 * self.axis.length

[docs]    def contains(self, point):
"""Test region membership of point in :class:Cylinder.

Parameters
----------
point : array_like

Returns
-------
:class:~python:bool
True if point is within :class:Cylinder,
False otherwise.

Notes
-----
A point :math:(p_x, p_y, p_z) is within the bounded region
of a cylinder with radius :math:r around the line from
:math:p_1=(x_1, y_1, z_1) to :math:p_2 = (x_2, y_2, z_2)
if the following is true:

.. math::

0\\le q\\le 1\\land
(x_1 - p_x + (x_2 - x_1) q)^2 +
(y_1 - p_y + (y_2 - y_1) q)^2 +
(z_1 - p_z + (z_2 - z_1) q)^2 \\le r^2

where :math:q is:

.. math::

q = \\frac{(p_x - x_1)(x_2 - x_1) +
(p_y - y_1)(y_2 - y_1) +
(p_z - z_1)(z_2 - z_1)}{(x_2 - x_1)^2 +
(y_2 - y_1)^2 + (z_2 - z_1)^2}

"""
px, py, pz = Point(point)
x1, y1, z1 = self.p1
x2, y2, z2 = self.p2
r = self.r

q1 = ((px - x1) * (x2 - x1) +
(py - y1) * (y2 - y1) +
(pz - z1) * (z2 - z1)) / \
((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2)

q2 = (x1 - px + (x2 - x1) * q1) ** 2 + \
(y1 - py + (y2 - y1) * q1) ** 2 + \
(z1 - pz + (z2 - z1) * q1) ** 2

return (not np.allclose(x1, x2) or not np.allclose(y1, y2) or
not np.allclose(z1, z2)) and r > 0 and q1 >= 0 and q1 <= 1 \
and q2 <= r ** 2

[docs]    def todict(self):
"""Returns a :class:~python:dict of the :class:Cylinder \
constructor parameters."""
return dict(p1=self.p1, p2=self.p2, r=self.r)

[docs]class Cone(Geometric3DRegion):
"""Geometric3DRegion for a cone.

Represents a cone with a base of radius :math:r centered at
:math:p_1=(x_1,y_1,z_1) and a tip at :math:p_2=(x_2, y_2, z_2).

Parameters
----------
p1, p2 : array_like, optional
3-tuples or :class:~sknano.core.Point class instances
for a :class:Cone with a base of radius r centered at
:math:p_1=(x_1,y_1,z_1) and a tip at :math:p_2=(x_2,y_2,z_2).
r : float, optional
Radius :math:r of :class:Cone base

Notes
-----
:class:Cone represents a bounded cone region
:math:\\left\\{p_1+\\rho(1-z)\\cos(\\theta)\\mathbf{v}_1 +
\\rho(1-z)\\sin(\\theta)\\mathbf{v}_2+\\mathbf{v}_3 z|
0\\le\\theta\\le 2\pi\\land
0\\le\\rho\\le 1\\land
0\\le z\\le 1\\right\\}
where :math:\\mathbf{v}_3=p_2-p_1 and vectors
:math:(\\mathbf{v}_1,\\mathbf{v}_2,\\mathbf{v}_3) are orthogonal
with :math:|\\mathbf{v}_1|=|\\mathbf{v}_2|=1 and
:math:p_1=(x_1,y_1,z_1) and :math:p_2=(x_2, y_2, z_2).

Calling :class:Cone with no parameters is equivalent to
:class:Cone\ (p1=[0, 0, 0], p2=[0, 0, 2], r=1).

"""
def __init__(self, p1=None, p2=None, r=1):
super().__init__()

self._p1 = Point()
self._p2 = Point()

if p1 is None:
p1 = [0, 0, 0]
if p2 is None:
p2 = [0, 0, 2]

self.p1 = p1
self.p2 = p2
self.points.extend([self.p1, self.p2])
self.r = r
self.fmtstr = "p1={p1!r}, p2={p2!r}, r={r:.3f}"

@property
def r(self):
"""Radius :math:r of :class:Cone base."""
return self._r

@r.setter
def r(self, value):
self._r = float(value)

@property
def p1(self):
"""Center point :math:(x_1, y_1, z_1) of :class:Cone base."""
return self._p1

@p1.setter
def p1(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._p1[:] = Point(value)

@property
def p2(self):
"""Point :math:(x_2, y_2, z_2) of :class:Cone tip."""
return self._p2

@p2.setter
def p2(self, value):
if not isinstance(value, (tuple, list, np.ndarray)) or len(value) != 3:
raise TypeError('Expected a 3-element array_like object')
self._p2[:] = Point(value)

@property
def axis(self):
""":class:Cone axis :class:Vector from \
:math:\\boldsymbol{\\ell}=p_2 - p_1

Returns
-------
:class:Vector

"""
return Vector(p0=self.p1, p=self.p2)

@property
def centroid(self):
""":class:Cone centroid, :math:(c_x, c_y, c_z).

Computed as:

.. math::

c_x = \\frac{3 x_1 + x_2}{4}

c_y = \\frac{3 y_1 + y_2}{4}

c_z = \\frac{3 z_1 + z_2}{4}

Returns
-------
:class:Point
3D :class:Point of :class:Cone centroid.

"""
return (3 * self.p1 + self.p2) / 4

@property
def volume(self):
""":class:Cone volume, :math:V=\\frac{1}{3}\\pi r^2 \\ell."""
return np.pi * self.r ** 2 * self.axis.length / 3

[docs]    def contains(self, point):
"""Test region membership of point in :class:Cone.

Parameters
----------
point : array_like

Returns
-------
:class:~python:bool
True if point is within :class:Cone,
False otherwise.

Notes
-----
A point :math:(p_x, p_y, p_z) is within the bounded region
of a cone with a base of radius :math:r centered at
:math:p_1=(x_1, y_1, z_1) and tip at :math:p_2 = (x_2, y_2, z_2)
if the following is true:

.. math::

0\\le q\\le 1\\land
(x_1 - p_x + (x_2 - x_1) q)^2 +
(y_1 - p_y + (y_2 - y_1) q)^2 +
(z_1 - p_z + (z_2 - z_1) q)^2 \\le r^2 q^2

where :math:q is:

.. math::

q = \\frac{(p_x - x_1)(x_2 - x_1) +
(p_y - y_1)(y_2 - y_1) +
(p_z - z_1)(z_2 - z_1)}{(x_2 - x_1)^2 +
(y_2 - y_1)^2 + (z_2 - z_1)^2}

"""
px, py, pz = Point(point)
x1, y1, z1 = self.p1
x2, y2, z2 = self.p2
r = self.r

q1 = ((px - x1) * (x2 - x1) +
(py - y1) * (y2 - y1) +
(pz - z1) * (z2 - z1)) / \
((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2)

q2 = (x1 - px + (x2 - x1) * q1) ** 2 + \
(y1 - py + (y2 - y1) * q1) ** 2 + \
(z1 - pz + (z2 - z1) * q1) ** 2

q3 = r ** 2 * q1 ** 2

return (not np.allclose(x1, x2) or not np.allclose(y1, y2) or
not np.allclose(z1, z2)) and r > 0 and q1 >= 0 and q1 <= 1 \
and q2 <= q3

[docs]    def todict(self):
"""Returns a :class:~python:dict of the :class:Cone \
constructor parameters."""
return dict(p1=self.p1, p2=self.p2, r=self.r)