Source code for skspatial.objects.plane

"""Module for the Plane class."""
from __future__ import annotations

from typing import Optional
from typing import Tuple

import numpy as np
from mpl_toolkits.mplot3d import Axes3D

from skspatial.objects._base_line_plane import _BaseLinePlane
from skspatial.objects._mixins import _ToPointsMixin
from skspatial.objects.line import Line
from skspatial.objects.point import Point
from skspatial.objects.points import Points
from skspatial.objects.vector import Vector
from skspatial.typing import array_like


[docs]class Plane(_BaseLinePlane, _ToPointsMixin): """ A plane in space. The plane is defined by a point and a normal vector. Parameters ---------- point : array_like Point on the plane. direction : array_like Normal vector of the plane. kwargs : dict, optional Additional keywords passed to :meth:`Vector.is_zero`. This method is used to ensure that the normal vector is not the zero vector. Attributes ---------- point : Point Point on the plane. normal : Vector Unit normal vector. vector : Vector Same as the normal. dimension : int Dimension of the plane. Raises ------ ValueError If the point and vector have different dimensions. If the vector is all zeros. Examples -------- >>> from skspatial.objects import Plane >>> plane = Plane(point=[0, 0, 0], normal=[0, 0, 5]) >>> plane Plane(point=Point([0, 0, 0]), normal=Vector([0, 0, 5])) >>> plane.normal Vector([0, 0, 5]) The normal can also be accessed with the ``vector`` attribute. >>> plane.vector Vector([0, 0, 5]) The plane dimension is the dimension of the point and vector. >>> plane.dimension 3 >>> Plane([0, 0], [1, 0, 0]) Traceback (most recent call last): ... ValueError: The point and vector must have the same dimension. >>> Plane([1, 1], [0, 0]) Traceback (most recent call last): ... ValueError: The vector must not be the zero vector. """ def __init__(self, point: array_like, normal: array_like): super().__init__(point, normal) self.normal = self.vector
[docs] @classmethod def from_vectors(cls, point: array_like, vector_a: array_like, vector_b: array_like, **kwargs) -> Plane: """ Instantiate a plane from a point and two vectors. The two vectors span the plane. Parameters ---------- point : array_like Point on the plane. vector_a, vector_b : array_like Input vectors. kwargs : dict, optional Additional keywords passed to :meth:`Vector.is_parallel`. Returns ------- Plane Plane containing input point and spanned by the two input vectors. Raises ------ ValueError If the vectors are parallel. Examples -------- >>> from skspatial.objects import Plane >>> Plane.from_vectors([0, 0], [1, 0], [0, 1]) Plane(point=Point([0, 0, 0]), normal=Vector([0, 0, 1])) >>> Plane.from_vectors([0, 0], [1, 0], [2, 0]) Traceback (most recent call last): ... ValueError: The vectors must not be parallel. """ vector_a = Vector(vector_a) if vector_a.is_parallel(vector_b, **kwargs): raise ValueError("The vectors must not be parallel.") # The cross product returns a 3D vector. vector_normal = vector_a.cross(vector_b) # Convert the point to 3D so that it matches the vector dimension. point = Point(point).set_dimension(3) return cls(point, vector_normal)
[docs] @classmethod def from_points(cls, point_a: array_like, point_b: array_like, point_c: array_like, **kwargs) -> Plane: """ Instantiate a plane from three points. The three points lie on the plane. Parameters ---------- point_a, point_b, point_c: array_like Three points defining the plane. kwargs: dict, optional Additional keywords passed to :meth:`Points.are_collinear`. Returns ------- Plane Plane containing the three input points. Raises ------ ValueError If the points are collinear. Examples -------- >>> from skspatial.objects import Plane >>> Plane.from_points([0, 0], [1, 0], [3, 3]) Plane(point=Point([0, 0, 0]), normal=Vector([0, 0, 3])) The order of the points affects the direction of the normal vector. >>> Plane.from_points([0, 0], [3, 3], [1, 0]) Plane(point=Point([0, 0, 0]), normal=Vector([ 0, 0, -3])) >>> Plane.from_points([0, 0], [0, 1], [0, 3]) Traceback (most recent call last): ... ValueError: The points must not be collinear. """ if Points([point_a, point_b, point_c]).are_collinear(**kwargs): raise ValueError("The points must not be collinear.") vector_ab = Vector.from_points(point_a, point_b) vector_ac = Vector.from_points(point_a, point_c) return Plane.from_vectors(point_a, vector_ab, vector_ac)
[docs] def cartesian(self) -> Tuple[np.number, np.number, np.number, np.number]: """ Return the coefficients of the Cartesian equation of the plane. The equation has the form ax + by + cz + d = 0. Returns ------- tuple Coefficients a, b, c, d. Raises ------ ValueError If the plane dimension is higher than 3. Examples -------- >>> from skspatial.objects import Plane >>> Plane(point=[1, 1], normal=[1, 0]).cartesian() (1, 0, 0, -1) >>> Plane(point=[1, 2, 0], normal=[0, 0, 1]).cartesian() (0, 0, 1, 0) >>> Plane(point=[1, 2, 8], normal=[0, 0, 5]).cartesian() (0, 0, 5, -40) >>> Plane(point=[4, 9, -1], normal=[10, 2, 4]).cartesian() (10, 2, 4, -54) >>> Plane([0, 0, 0, 0], [1, 0, 0, 0]).cartesian() Traceback (most recent call last): ... ValueError: The plane dimension must be <= 3. """ if self.dimension > 3: raise ValueError("The plane dimension must be <= 3.") # The normal must be 3D to extract the coefficients. a, b, c = self.normal.set_dimension(3) d = -self.normal.dot(self.point) return a, b, c, d
[docs] def project_point(self, point: array_like) -> Point: """ Project a point onto the plane. Parameters ---------- point : array_like Input point. Returns ------- Point Projection of the point onto the plane. Examples -------- >>> from skspatial.objects import Plane >>> plane = Plane(point=[0, 0, 0], normal=[0, 0, 2]) >>> plane.project_point([10, 2, 5]) Point([10., 2., 0.]) >>> plane.project_point([5, 9, -3]) Point([5., 9., 0.]) """ # Vector from the point in space to the point on the plane. vector_to_plane = Vector.from_points(point, self.point) # Perpendicular vector from the point in space to the plane. vector_projected = self.normal.project_vector(vector_to_plane) return Point(point) + vector_projected
[docs] def project_vector(self, vector: array_like) -> Vector: """ Project a vector onto the plane. Parameters ---------- vector : array_like Input vector. Returns ------- Vector Projection of the vector onto the plane. Examples -------- >>> from skspatial.objects import Plane >>> plane = Plane([0, 4, 0], [0, 1, 1]) >>> plane.project_vector([2, 4, 8]) Vector([ 2., -2., 2.]) """ point_in_space = self.point + vector point_on_plane = self.project_point(point_in_space) return Vector.from_points(self.point, point_on_plane)
[docs] def project_line(self, line: Line, **kwargs: float) -> Line: """ Project a line onto the plane. This method can also handle the case where the line is parallel to the plane. Parameters ---------- line : Line Input line. kwargs : dict, optional Additional keywords passed to :meth:`Vector.is_perpendicular`, which is used to check if the line is parallel to the plane (i.e., the line direction is perpendicular to the plane normal). Returns ------- Line Projection of the line onto the plane. Examples -------- >>> from skspatial.objects import Line, Plane >>> plane = Plane([0, 0, 0], [0, 0, 1]) >>> line = Line([0, 0, 0], [1, 1, 1]) >>> plane.project_line(line) Line(point=Point([0., 0., 0.]), direction=Vector([1., 1., 0.])) The line is parallel to the plane. >>> line = Line([0, 0, 5], [1, 0, 0]) >>> plane.project_line(line) Line(point=Point([0., 0., 0.]), direction=Vector([1, 0, 0])) """ if self.normal.is_parallel(line.vector, **kwargs): raise ValueError("The line and plane must not be perpendicular.") point_projected = self.project_point(line.point) if self.normal.is_perpendicular(line.vector, **kwargs): return Line(point_projected, line.vector) vector_projected = self.project_vector(line.vector) return Line(point_projected, vector_projected)
[docs] def distance_point_signed(self, point: array_like) -> np.float64: """ Return the signed distance from a point to the plane. Parameters ---------- point : array_like Input point. Returns ------- np.float64 Signed distance from the point to the plane. References ---------- http://mathworld.wolfram.com/Point-PlaneDistance.html Examples -------- >>> from skspatial.objects import Plane >>> plane = Plane([0, 0, 0], [0, 0, 1]) >>> plane.distance_point_signed([5, 2, 0]) 0.0 >>> plane.distance_point_signed([5, 2, 1]) 1.0 >>> plane.distance_point([5, 2, -4]) 4.0 >>> plane.distance_point_signed([5, 2, -4]) -4.0 """ vector_to_point = Vector.from_points(self.point, point) return self.normal.scalar_projection(vector_to_point)
def distance_point(self, point: array_like) -> np.float64: """ Return the distance from a point to the plane. Parameters ---------- point : array_like Input point. Returns ------- np.float64 Distance from the point to the plane. References ---------- http://mathworld.wolfram.com/Point-PlaneDistance.html Examples -------- >>> from skspatial.objects import Plane >>> plane = Plane([0, 0, 0], [0, 0, 1]) >>> plane.distance_point([5, 2, 0]) 0.0 >>> plane.distance_point([5, 2, 1]) 1.0 >>> plane.distance_point([5, 2, -4]) 4.0 """ return abs(self.distance_point_signed(point))
[docs] def side_point(self, point: array_like) -> int: """ Find the side of the plane where a point lies. Parameters ---------- point : array_like Input point. Returns ------- int -1 if the point is behind the plane. 0 if the point is on the plane. 1 if the point is in front of the plane. Examples -------- >>> from skspatial.objects import Plane >>> plane = Plane([0, 0, 0], [0, 0, 1]) The point is in on the plane. >>> plane.side_point([2, 5, 0]) 0 The point is in front of the plane. >>> plane.side_point([1, -5, 6]) 1 The point is behind the plane. >>> plane.side_point([5, 8, -4]) -1 Higher dimensions are supported. >>> plane = Plane([0, 0, 0, 0], [0, 1, 0, 1]) >>> plane.side_point([0, -10, 4, 1]) -1 """ return int(np.sign(self.distance_point_signed(point)))
[docs] def intersect_line(self, line: Line, **kwargs) -> Point: """ Intersect the plane with a line. The line and plane must not be parallel. Parameters ---------- line : Line Input line. kwargs : dict, optional Additional keywords passed to :meth:`Vector.is_perpendicular`. Returns ------- Point The point of intersection. Raises ------ ValueError If the line and plane are parallel. References ---------- http://geomalgorithms.com/a05-_intersect-1.html Examples -------- >>> from skspatial.objects import Line, Plane >>> line = Line([0, 0, 0], [0, 0, 1]) >>> plane = Plane([0, 0, 0], [0, 0, 1]) >>> plane.intersect_line(line) Point([0., 0., 0.]) >>> plane = Plane([2, -53, -7], [0, 0, 1]) >>> plane.intersect_line(line) Point([ 0., 0., -7.]) >>> line = Line([0, 1, 0], [1, 0, 0]) >>> plane.intersect_line(line) Traceback (most recent call last): ... ValueError: The line and plane must not be parallel. """ if self.normal.is_perpendicular(line.direction, **kwargs): raise ValueError("The line and plane must not be parallel.") vector_plane_line = Vector.from_points(self.point, line.point) num = -self.normal.dot(vector_plane_line) denom = self.normal.dot(line.direction) # Vector along the line to the intersection point. vector_line_scaled = num / denom * line.direction return line.point + vector_line_scaled
[docs] def intersect_plane(self, other: Plane, **kwargs) -> Line: """ Intersect the plane with another. The planes must not be parallel. Parameters ---------- other : Plane Other plane. kwargs : dict, optional Additional keywords passed to :meth:`Vector.is_parallel`. Returns ------- Line The line of intersection. Raises ------ ValueError If the planes are parallel. References ---------- http://tbirdal.blogspot.com/2016/10/a-better-approach-to-plane-intersection.html Examples -------- >>> from skspatial.objects import Plane >>> plane_a = Plane([0, 0, 0], [0, 0, 1]) >>> plane_b = Plane([0, 0, 0], [1, 0, 0]) >>> plane_a.intersect_plane(plane_b) Line(point=Point([0., 0., 0.]), direction=Vector([0, 1, 0])) >>> plane_b = Plane([5, 16, -94], [1, 0, 0]) >>> plane_a.intersect_plane(plane_b) Line(point=Point([5., 0., 0.]), direction=Vector([0, 1, 0])) >>> plane_b = Plane([0, 0, 1], [1, 0, 1]) >>> plane_a.intersect_plane(plane_b) Line(point=Point([1., 0., 0.]), direction=Vector([0, 1, 0])) >>> plane_b = Plane([0, 0, 5], [0, 0, -8]) >>> plane_a.intersect_plane(plane_b) Traceback (most recent call last): ... ValueError: The planes must not be parallel. """ if self.normal.is_parallel(other.normal, **kwargs): raise ValueError("The planes must not be parallel.") array_normals_stacked = np.vstack((self.normal, other.normal)) # Construct a matrix for a linear system. array_00 = 2 * np.eye(3) array_01 = array_normals_stacked.T array_10 = array_normals_stacked array_11 = np.zeros((2, 2)) matrix = np.block([[array_00, array_01], [array_10, array_11]]) dot_a = np.dot(self.point, self.normal) dot_b = np.dot(other.point, other.normal) array_y = np.array([0, 0, 0, dot_a, dot_b]) # Solve the linear system. solution = np.linalg.solve(matrix, array_y) point_line = Point(solution[:3]) direction_line = self.normal.cross(other.normal) return Line(point_line, direction_line)
[docs] @classmethod def best_fit(cls, points: array_like, tol: Optional[float] = None, **kwargs) -> Plane: """ Return the plane of best fit for a set of 3D points. Parameters ---------- points : array_like Input 3D points. tol : float | None, optional Keyword passed to :meth:`Points.are_collinear` (default None). kwargs : dict, optional Additional keywords passed to :func:`numpy.linalg.svd` Returns ------- Plane The plane of best fit. Raises ------ ValueError If the points are collinear or are not 3D. References ---------- Using SVD for some fitting problems Inge Söderkvist Algorithm 3.1 https://www.ltu.se/cms_fs/1.51590!/svd-fitting.pdf Examples -------- >>> from skspatial.objects import Plane >>> points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]] >>> plane = Plane.best_fit(points) The point on the plane is the centroid of the points. >>> plane.point Point([0.25, 0.25, 0.25]) The plane normal is a unit vector. >>> plane.normal.round(3) Vector([-0.577, -0.577, -0.577]) >>> points = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0]] >>> Plane.best_fit(points) Plane(point=Point([0.5, 0.5, 0. ]), normal=Vector([0., 0., 1.])) >>> Plane.best_fit(points, full_matrices=False) Plane(point=Point([0.5, 0.5, 0. ]), normal=Vector([0., 0., 1.])) """ points = Points(points) if points.dimension != 3: raise ValueError("The points must be 3D.") if points.are_collinear(tol=tol): raise ValueError("The points must not be collinear.") points_centered, centroid = points.mean_center(return_centroid=True) u, _, _ = np.linalg.svd(points_centered.T, **kwargs) normal = Vector(u[:, 2]) return cls(centroid, normal)
[docs] def to_mesh( self, lims_x: array_like = (-1, 1), lims_y: array_like = (-1, 1), ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ Return coordinate matrices for the 3D surface of the plane. Parameters ---------- lims_x, lims_y : (2,) tuple x and y limits of the plane. Tuple of form (min, max). The default is (-1, 1). The point on the plane is used as the origin. Returns ------- X, Y, Z: ndarray Coordinate matrices. Examples -------- >>> from skspatial.objects import Plane >>> X, Y, Z = Plane([0, 0, 0], [0, 0, 1]).to_mesh() >>> X array([[-1, 1], [-1, 1]]) >>> Y array([[-1, -1], [ 1, 1]]) >>> Z array([[0., 0.], [0., 0.]]) """ a, b, c, d = self.cartesian() x_center, y_center = self.point[:2] values_x = x_center + lims_x values_y = y_center + lims_y X, Y = np.meshgrid(values_x, values_y) if c != 0: Z = -(a * X + b * Y + d) / c elif b != 0: Z = -(a * X + c * Y + d) / b X, Y, Z = X, Z, Y else: Z = -(b * X + c * Y + d) / a X, Y, Z = Z, X, Y return X, Y, Z
[docs] def plot_3d(self, ax_3d: Axes3D, lims_x: array_like = (-1, 1), lims_y: array_like = (-1, 1), **kwargs) -> None: """ Plot a 3D plane. Parameters ---------- ax_3d : Axes3D Instance of :class:`~mpl_toolkits.mplot3d.axes3d.Axes3D`. lims_x, lims_y : (2,) tuple x and y limits of the plane. Tuple of form (min, max). The default is (-1, 1). The point on the plane is used as the origin. kwargs : dict, optional Additional keywords passed to :meth:`~mpl_toolkits.mplot3d.axes3d.Axes3D.plot_surface`. Examples -------- .. plot:: :include-source: >>> import matplotlib.pyplot as plt >>> from mpl_toolkits.mplot3d import Axes3D >>> from skspatial.objects import Plane >>> fig = plt.figure() >>> ax = fig.add_subplot(111, projection='3d') >>> plane = Plane([5, 3, 1], [1, 0, 1]) >>> plane.plot_3d(ax, alpha=0.2) >>> plane.point.plot_3d(ax, s=100) """ X, Y, Z = self.to_mesh(lims_x, lims_y) ax_3d.plot_surface(X, Y, Z, **kwargs)