"""Module for converting visualizable objects to discrete array based meshes or curves.
"""
import dataclasses
from collections.abc import Sequence
from functools import partial
import numpy as np
import ddg
import ddg._error_messages as em
_default_curve_sampling = (100, 0.1)
_default_surface_sampling = (30, 0.1)
_default_subspace_size = 100
_default_curve_radius = 0.015
_default_point_radius = 0.05
_default_icosphere_subdivision_steps = 3
def _default_point_mesh(
subdivision_steps=_default_icosphere_subdivision_steps, radius=_default_point_radius
):
"""Create default point mesh"""
return from_half_edge_surface(
ddg.halfedge.icosphere(
subdivision_steps,
radius,
)
)
def _from_same_dimensional_components(components, /):
"""Combine same dimensional components.
Parameters
----------
components : Sequence[numpy.ndarray, Curve, Mesh]
Returns
-------
Points | CurveList | Mesh | components
Returns Points if all components are instances Points.
Returns DisconnectedCurve if all components are instances of Curve.
Returns Mesh if all components are instances of Mesh.
Returns components otherwise.
"""
if len(components) == 0:
return Mesh(np.zeros((0, 3), float), ())
if all(isinstance(c, Points) for c in components):
return Points(np.vstack([c.points for c in components]))
elif all(isinstance(c, Curve) for c in components):
return CurveList(components)
elif all(isinstance(c, Mesh) for c in components):
return mesh_union(components)
else:
return components
[docs]def from_discrete_net(net, /):
"""Compute all points and faces of a discrete net.
Parameters
----------
net : ddg.nets.DiscreteNet
A net with two-dimensional domain and `net.domain.bounded == True`.
Returns
-------
Mesh
See Also
--------
ddg.arrays.from_discrete_curve
Examples
--------
>>> import ddg
>>> net = ddg.nets.DiscreteNet(lambda u, v: (u, v), ((0, 2), (-4, -1)))
>>> points, faces, _ = ddg.arrays.from_discrete_net(net)
>>> points
array([[ 0, -4],
[ 0, -3],
[ 0, -2],
[ 0, -1],
[ 1, -4],
[ 1, -3],
[ 1, -2],
[ 1, -1],
[ 2, -4],
[ 2, -3],
[ 2, -2],
[ 2, -1]])
>>> faces
array([[ 0, 4, 5, 1],
[ 1, 5, 6, 2],
[ 2, 6, 7, 3],
[ 4, 8, 9, 5],
[ 5, 9, 10, 6],
[ 6, 10, 11, 7]])
"""
if not net.domain.bounded:
raise ValueError(
"net.domain.bounded must be True but "
f"net.domain.bounded = {net.domain.bounded}"
)
return Mesh(np.array(list(net[net.domain])), np.array(list(net.domain.face_data)))
[docs]def from_discrete_curve(curve, /):
"""Compute all points of a discrete curve.
Parameters
----------
curve : ddg.nets.DiscreteCurve
A curve with `curve.domain.bounded == True`.
Returns
-------
Curve
See Also
--------
ddg.arrays.from_discrete_net
Examples
--------
>>> import ddg
>>> curve = ddg.nets.DiscreteCurve(lambda t: (t,), (-1, 1, False))
>>> points, periodicity = ddg.arrays.from_discrete_curve(curve)
>>> points
array([[-1],
[ 0],
[ 1]])
>>> periodicity
False
"""
if not curve.domain.bounded:
raise ValueError(
f"curve.domain.bounded must be True but {curve.domain.bounded = }"
)
return Curve(np.array(list(curve[curve.domain])), curve.domain.periodic)
[docs]def from_empty_net(_, /):
"""Convert an empty net to an empty mesh.
Parameters
----------
_ : ddg.nets.EmptyNet
Returns
-------
Mesh
Examples
--------
>>> import ddg
>>> empty_net = ddg.nets.EmptyNet(())
>>> points, faces, _ = ddg.arrays.from_empty_net(empty_net)
>>> points
array([], shape=(0, 3), dtype=float64)
>>> faces
()
"""
return Mesh(np.zeros((0, 3), float), ())
[docs]def from_point_net(net, /):
"""Convert a point net to a mesh of one point and no faces.
Parameters
----------
net : ddg.nets.PointNet
Returns
-------
Points
The object contains the value of `net()`, which is a point in R^n.
Examples
--------
>>> import ddg
>>> points_net = ddg.nets.PointNet((1, 2, 3))
>>> ddg.arrays.from_point_net(points_net)
Points(points=array([[1, 2, 3]]))
"""
return Points(np.array(net()))
def _from_net(
net,
curve_sampling=_default_curve_sampling,
surface_sampling=_default_surface_sampling,
):
if isinstance(net, ddg.nets.DiscreteCurve):
return from_discrete_curve(net)
elif isinstance(net, ddg.nets.DiscreteNet):
return from_discrete_net(net)
elif isinstance(net, ddg.nets.EmptyNet):
return from_empty_net(net)
elif isinstance(net, ddg.nets.PointNet):
return from_point_net(net)
elif isinstance(net, ddg.nets.SmoothCurve):
sampling = [curve_sampling[1], curve_sampling[0], "c"]
dnet = ddg.nets.sample_smooth_net(net, sampling)
return _from_net(dnet)
elif isinstance(net, ddg.nets.SmoothNet):
sampling = [surface_sampling[1], surface_sampling[0], "c"]
dnet = ddg.nets.sample_smooth_net(net, sampling)
return _from_net(dnet)
else:
raise Exception(em.should_never_happen)
[docs]def from_net_collection(
net_collection,
/,
curve_sampling=_default_curve_sampling,
surface_sampling=_default_surface_sampling,
):
"""Convert a net collection to a (disconnected) curve or mesh.
Parameters
----------
nets : ddg.nets.NetCollection
Every net in `nets` has to be an instance of the same of one of these::
- a SmoothCurve
- a SmoothNet
- a DiscreteCurve with bounded domain
- a DiscreteNet with bounded domain
- a PointNet
- an EmptyNet
Additionally, the dimensions of the images each net in `nets` must be
equal to one another.
curve_sampling : tuple of an int and a float
(default=_default_curve_sampling)
Determines the sample number and stepsize for curves
surface_sampling : tuple of an int and a float
(default=_default_surface_sampling)
Determines the sample number and stepsize for surfaces
Returns
-------
Points, CurveList or Mesh
Raises
------
ValueError
If the nets contained in `nets` don't satisfy the conditions in the
description of `nets`.
Examples
--------
>>> import ddg
>>> from ddg.arrays import from_net_collection
>>> curve_0 = ddg.nets.DiscreteCurve(lambda t: (t, 0), (0, 2, False))
>>> curve_1 = ddg.nets.DiscreteCurve(lambda t: (t, t**2), (0, 2, True))
>>> net_collection = ddg.nets.NetCollection([curve_0, curve_1])
>>> (points_0, periodicity_0), (points_1, periodicity_1) = from_net_collection(
... net_collection
... )
>>> points_0
array([[0, 0],
[1, 0],
[2, 0]])
>>> periodicity_0
False
>>> points_1
array([[0, 0],
[1, 1],
[2, 4]])
>>> periodicity_1
True
>>> two_sheeted_hyperboloid = ddg.geometry.Quadric(np.diag([1, 1, -1, 1]))
>>> smooth_net_collection = ddg.to_smooth_net(two_sheeted_hyperboloid)
>>> mesh = ddg.arrays.from_net_collection(
... smooth_net_collection, surface_sampling=(3, 0.1)
... )
>>> mesh.points
array([[ 0. , 0. , -1. ],
[ 0. , 0. , -1. ],
[ 0. , 0. , -1. ],
[ 0. , -0.10016675, -1.00500417],
[ 0.08674695, 0.05008338, -1.00500417],
[-0.08674695, 0.05008338, -1.00500417],
[ 0. , -0.201336 , -1.02006676],
[ 0.17436209, 0.100668 , -1.02006676],
[-0.17436209, 0.100668 , -1.02006676],
[ 0. , 0. , 1. ],
[ 0. , 0. , 1. ],
[ 0. , 0. , 1. ],
[ 0. , -0.10016675, 1.00500417],
[ 0.08674695, 0.05008338, 1.00500417],
[-0.08674695, 0.05008338, 1.00500417],
[ 0. , -0.201336 , 1.02006676],
[ 0.17436209, 0.100668 , 1.02006676],
[-0.17436209, 0.100668 , 1.02006676]])
>>> mesh.faces
array([[ 0, 3, 4, 1],
[ 1, 4, 5, 2],
[ 3, 6, 7, 4],
[ 4, 7, 8, 5],
[ 0, 2, 5, 3],
[ 3, 5, 8, 6],
[ 9, 12, 13, 10],
[10, 13, 14, 11],
[12, 15, 16, 13],
[13, 16, 17, 14],
[ 9, 11, 14, 12],
[12, 14, 17, 15]])
"""
try:
net_collection.dimension
except AttributeError as ae:
if (
"Dimension of NetCollection is not well-defined because it "
"contains nets of dimensions " in str(ae)
):
raise ValueError(
"The NetCollection contains nets of different dimensions. "
"Consider converting the nets separately."
) from ae
raise AttributeError(ae)
# NOTE: This also catches the case where an empty net collection is given.
# It will then be converted to an empty mesh
if isinstance(net_collection[0], ddg.nets.SmoothNet):
if net_collection.dimension == 1:
sampling = [curve_sampling[1], curve_sampling[0], "c"]
else:
sampling = [surface_sampling[1], surface_sampling[0], "c"]
net_collection = ddg.nets.sample_smooth_net(net_collection, sampling)
def convertible(net):
if isinstance(net, ddg.nets.DiscreteNet):
if net.domain.bounded:
return net
else:
return "The net is a DiscreteNet, but its domain isn't bounded"
elif isinstance(net, (ddg.nets.EmptyNet, ddg.nets.PointNet)):
return net
else:
return (
f"The net's type is {type(net)}, so it is neither a DiscreteCurve "
"or DiscreteNet with a bounded domain, an EmptyNet or a PointNet"
)
# It is very important to use ._nets rather than .nets because net
# collections maintain two lists of nets for some reason
nets_or_errors = list(map(convertible, net_collection._nets))
bad = [
(i, error) for i, error in enumerate(nets_or_errors) if isinstance(error, str)
]
if bad:
raise ValueError(
"Some nets in the collection can't be converted to arrays:\n"
+ "\n".join(f"Net {i}: {error}" for i, error in bad)
)
return _from_same_dimensional_components(list(map(_from_net, nets_or_errors)))
def _he_verts(surface, coordinate_attribute):
return np.array([getattr(v, coordinate_attribute) for v in surface.verts])
def _isolated_he_edges(surface):
edges = [
e
for e in ddg.halfedge.single_edges(surface)
if e.face is None and e.opp.face is None
]
if edges:
return np.array([(e.head.index, e.pre.head.index) for e in edges])
else:
return np.zeros((0, 2), int)
def _he_faces(surface):
return [
tuple(v.index for v in ddg.halfedge.face_vertices(f)) for f in surface.faces
]
[docs]def from_half_edge_surface(surface, coordinate_attribute="co", /):
"""Get the points, isolated edges and faces of a half edge surface.
Parameters
----------
surface : ddg.halfedge.Surface
coordinate_attribute : str (default="co")
Returns
-------
Mesh
Examples
--------
>>> import ddg
>>> square = ddg.halfedge.grid((2, 2))
>>> (
... points,
... faces,
... isolated_edges_square,
... ) = ddg.arrays.from_half_edge_surface(square)
>>> points
array([[0, 0],
[1, 0],
[0, 1],
[1, 1]])
>>> faces
array([[1, 3, 2, 0]])
>>> isolated_edges_square
array([], shape=(0, 2), dtype=int64)
>>> surface = ddg.halfedge.grid((2, 2))
>>> vertex = ddg.halfedge.some_vertex(surface)
>>> isolated_vertex = surface.add_vertex()
>>> isolated_vertex.co = np.array([2, 2])
>>> surface.add_edge(vertex, isolated_vertex)
Edges(index=8, head=4, face=None, surface=...)
>>> (_, _, isolated_edges_surface) = ddg.arrays.from_half_edge_surface(surface)
>>> isolated_edges_surface
array([[4, 0]])
"""
return Mesh(
_he_verts(surface, coordinate_attribute),
_he_faces(surface),
_isolated_he_edges(surface),
)
[docs]def from_0d_subspace(point, /):
"""Convert a point subspace to a mesh of one point and no faces.
Parameters
----------
point : ddg.geometry.Subspace
A subspace with `point.dimension == 0`.
Returns
-------
Points
Examples
--------
>>> import ddg
>>> point = ddg.geometry.subspace_from_affine_points((0, 0))
>>> ddg.arrays.from_0d_subspace(point)
Points(points=array([[0., 0.]]))
"""
if point.dimension != 0:
message = em.build_error_message(
point,
"point",
None,
[em.str_(f"point.dimension must be 0 but {point.dimension = }")],
)
raise ValueError(message)
return Points(point.affine_point)
[docs]def from_1d_subspace(line, length, /):
"""Sample a 1-dimensional subspace.
Parameters
----------
line : ddg.geometry.Subspace
A subspace with `line.dimension == 1`.
length : float
Must be greater than or equal to 0.0.
Returns
-------
Curve(ab, False)
`ab` is an array_like of shape `(subspace.ambient_dimension, 2)`
Let `a = ab[0]` and `b = ab[1]`, then `b - a` has length `length` and
`b - a / 2` is the closest point of the line to the origin.
Raises
------
ValueError
If `line.dimension != 1` or `length < 0.0`.
Examples
--------
>>> import ddg
>>> line = ddg.geometry.subspace_from_affine_points((0, 0, 0), (1, 0, 0))
>>> ab, periodicity = ddg.arrays.from_1d_subspace(line, 10)
>>> ab
array([[-5., 0., 0.],
[ 5., 0., 0.]])
>>> periodicity
False
"""
if line.dimension != 1:
message = em.build_error_message(
line,
"line",
None,
[em.str_(f"line.dimension must be 1 but {line.dimension = }")],
)
raise ValueError(message)
if length < 0:
raise ValueError(f"{length = } must be non-negative.")
match line.affine_point_and_directions:
case p, [d]:
d_ = length / 2 * d
return Curve(np.array((p - d_, p + d_)), False)
case _:
raise Exception(em.should_never_happen)
[docs]def from_2d_subspace(plane, side_length_1, side_length_2, /):
"""Convert a 2-dimensional subspace to a mesh.
The mesh is given as two triangles::
d c
┌─────┐ ┐
│ /│ │
│ / │ │
│ x │ │ side_length_2
│ / │ │
│/ │ │
└─────┘ ┘
a b
└─────┘
side_length_1
such that `x` is the closest point of the subspace to the origin.
Parameters
----------
plane : ddg.geometry.Subspace
A subspace with `subspace.dimension == 2`.
side_length_1 : float
Must be greater than or equal to 0.0.
side_length_2 : float
Must be greater than or equal to 0.0.
Returns
-------
Mesh(points, faces)
The first entry points has shape (4, 3) and comprises of a, b, c, d.
The second entry faces comprises of two tuples of 3 integers encoding
the two triangles.
Raises
------
ValueError
If `plane.dimension != 2` or `side_length_1 < 0.0` or `side_length_2 < 0.0`.
Examples
--------
>>> import ddg
>>> plane = ddg.geometry.subspace_from_affine_points(
... (0, 0, 0), (1, 0, 0), (0, 1, 0)
... )
>>> points, faces, _ = ddg.arrays.from_2d_subspace(plane, 10, 4)
>>> points
array([[-5., -2., 0.],
[ 5., -2., 0.],
[ 5., 2., 0.],
[-5., 2., 0.]])
>>> faces
array([[0, 1, 2],
[0, 2, 3]])
"""
if plane.dimension != 2:
raise ValueError(f"{plane.dimension = } must be 2.")
if side_length_1 < 0 or side_length_2 < 0:
raise ValueError(
f"{side_length_1 = } and {side_length_2 = } must be non-negative."
)
match plane.affine_point_and_directions:
case p, [u, v]:
mu = side_length_1 / 2 * u
mv = side_length_2 / 2 * v
a = p - mu - mv
b = p + mu - mv
c = p + mu + mv
d = p - mu + mv
points = np.row_stack((a, b, c, d))
faces = (0, 1, 2), (0, 2, 3)
return Mesh(points, faces)
case _:
raise Exception(em.should_never_happen)
[docs]def from_quadric(quadric, /, *num_step):
"""Convert a quadric to a curve or a mesh or a list thereof.
Every convertible quadric has a parametrization, which is sampled. These
parametrizations are implementation details, but the following guarantees
hold:
- every parametrization is defined on a box I_1 x ... x I_d, where d is the
dimension of the quadric and I_1, ..., I_d are possibly unbounded
intervals
- the k-th `num_step` argument determines how I_k is sampled
- the first entry `num` of the k-th `num_step` determines the number of
samples of I_k
- if I_k is unbounded, the second entry `step` of the k-th `num_step`
determines the step size between the samples of I_k
Consequently,
- increasing the `num` arguments increases the number of points and
additionally the number of faces if the quadric is 2-dimensional.
- if I_k isn't bounded, decreasing the respective `step` results in a finer
curve or mesh.
Parameters
----------
quadric : ddg.geometry.Quadric
*num_step : tuple of an int and a float
The length of `num_step` must be 0 if `quadric.dimension` is -1 or 0
and otherwise must be equal to `quadric.dimension`.
Returns
-------
Points, Curve, CurveList or Mesh
Raises
------
ValueError
If the quadric is a (possibly empty) set of points and `num_step` is
given.
Examples
--------
>>> import numpy as np
>>> import ddg
>>> xaxis = ddg.geometry.Quadric(np.diag([0, 1, 0]))
>>> xaxis_points, xaxis_periodicity = ddg.arrays.from_quadric(xaxis, (4, 0.1))
>>> xaxis_points
array([[ 0.1, 0. ],
[ 0. , 0. ],
[-0.1, 0. ],
[-0.2, 0. ]])
>>> xaxis_periodicity
False
A coarser sampling
>>> xaxis_points_coarser, _ = ddg.arrays.from_quadric(xaxis, (4, 0.5))
>>> xaxis_points_coarser
array([[ 0.5, 0. ],
[ 0. , 0. ],
[-0.5, 0. ],
[-1. , 0. ]])
>>> two_sheeted_hyperboloid = ddg.geometry.Quadric(np.diag([1, 1, -1, 1]))
>>> mesh = ddg.arrays.from_quadric(two_sheeted_hyperboloid, (3, 0.1), (3, 0.1))
>>> mesh.points
array([[ 0. , 0. , -1. ],
[ 0. , 0. , -1. ],
[ 0. , 0. , -1. ],
[ 0. , -0.10016675, -1.00500417],
[ 0.08674695, 0.05008338, -1.00500417],
[-0.08674695, 0.05008338, -1.00500417],
[ 0. , -0.201336 , -1.02006676],
[ 0.17436209, 0.100668 , -1.02006676],
[-0.17436209, 0.100668 , -1.02006676],
[ 0. , 0. , 1. ],
[ 0. , 0. , 1. ],
[ 0. , 0. , 1. ],
[ 0. , -0.10016675, 1.00500417],
[ 0.08674695, 0.05008338, 1.00500417],
[-0.08674695, 0.05008338, 1.00500417],
[ 0. , -0.201336 , 1.02006676],
[ 0.17436209, 0.100668 , 1.02006676],
[-0.17436209, 0.100668 , 1.02006676]])
>>> mesh.faces
array([[ 0, 3, 4, 1],
[ 1, 4, 5, 2],
[ 3, 6, 7, 4],
[ 4, 7, 8, 5],
[ 0, 2, 5, 3],
[ 3, 5, 8, 6],
[ 9, 12, 13, 10],
[10, 13, 14, 11],
[12, 15, 16, 13],
[13, 16, 17, 14],
[ 9, 11, 14, 12],
[12, 14, 17, 15]])
"""
# We intentionally don't allow for other ways to sample the
# parametrization, because the choice of parametrization is
# arbitrary and therefore an implementation detail. For example, it
# doesn't make sense to let the user specify samples in the domain
# such as anchor values.
#
# The parametrization isn't even always obvious, e.g. the paraboloid
# may be parametrized by
# (u, v) -> (u, v, u^2 + v^2)
# and also as a surface of revolution and a circle can be
# parametrized in two directions. So if the user wants more control,
# they should just directly sample a parametrization of their
# choice.
# TODO: Rewrite this using vectorised parametrization and domains
# once they are available.
net_or_collection = ddg.conversion.nets.geometry.quadrics.quadric_to_smooth_net(
quadric
)
if quadric.dimension == -1:
if not num_step:
sampling = ()
else:
raise ValueError(
em.build_error_message(
quadric,
"quadric",
"The quadric is empty, so num_step must be empty.",
[em.str_(f"{quadric.dimension = }")],
)
+ em.build_error_message(num_step, "num_step", "", [em.len_])
)
elif quadric.dimension >= 0:
if quadric.dimension == len(num_step):
if quadric.dimension == 1:
((num, step),) = num_step
# sampling is not allowed to be [(step, num, "c")], otherwise
# sample_interval raises a ValueError for some reason.
# Otherwise, this if branch could be merged with its else branch.
sampling = (step, num, "c")
else:
sampling = [(step, num, "c") for num, step in num_step]
else:
text = (
"To sample a quadric of quadric.dimension >= 0, num_step must have\n"
"length quadric.dimension."
)
raise ValueError(
em.build_error_message(
quadric,
"quadric",
text,
[em.str_(f"{quadric.dimension = }")],
)
+ em.build_error_message(num_step, "num_step", "", [em.len_])
)
else:
raise ValueError(em.should_never_happen)
if isinstance(net_or_collection, ddg.nets.SmoothNet):
discrete_net = ddg.nets.sample_smooth_net(net_or_collection, sampling)
return _from_net(discrete_net)
elif isinstance(net_or_collection, (ddg.nets.EmptyNet, ddg.nets.PointNet)):
return _from_net(net_or_collection)
elif isinstance(net_or_collection, ddg.nets.NetCollection):
discrete_net_collection = ddg.nets.sample_smooth_net(
net_or_collection, sampling
)
return from_net_collection(discrete_net_collection)
else:
raise Exception(em.should_never_happen)
[docs]def from_quadric_sphere(sphere, num, /):
"""Convert a sphere to a curve or a mesh.
Circles are converted to curves and 2-dimensional spheres are
converted to meshes.
Parameters
----------
sphere : ddg.geometry.spheres.QuadricSphere
A sphere that can be represented as a quadric.
num : int
Increasing this number increases the number of samples in each
direction of the sphere parametrization.
Returns
-------
Curve or Mesh
A curve is an tuple of points of type array_like and the
periodicity of type bool.
A mesh is an tuple of points of type array_like and the
faces of type array_like.
Examples
--------
>>> import numpy as np
>>> import ddg
>>> euc = ddg.geometry.euclidean(3)
>>> s1 = euc.sphere((0, 0, 1), 1)
>>> s1_points, periodicity = ddg.arrays.from_quadric_sphere(s1, 4)
>>> s1_points
array([[-1.0000000e+00, 0.0000000e+00],
[-6.1232340e-17, -1.0000000e+00],
[ 1.0000000e+00, -1.2246468e-16],
[ 1.8369702e-16, 1.0000000e+00]])
>>> periodicity
True
>>> s2 = euc.sphere((0, 0, 0, 1), 1)
>>> s2_points, faces, _ = ddg.arrays.from_quadric_sphere(s2, 4)
>>> s2_points
array([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
...
[-1.00000000e+00, -7.49879891e-33, -1.22464680e-16],
[-1.00000000e+00, 1.22464680e-16, -1.49975978e-32],
[-1.00000000e+00, 2.24963967e-32, 1.22464680e-16]])
>>> faces
array([[ 0, 4, 5, 1],
[ 1, 5, 6, 2],
[ 2, 6, 7, 3],
...
[ 0, 3, 7, 4],
[ 4, 7, 11, 8],
[ 8, 11, 15, 12]])
"""
sampling = sphere.dimension * [(num, 0.1)]
return from_quadric(sphere.quadric(), *sampling)
[docs]def from_indexed_face_set(ifs, coordinate_key="co", /):
"""Get the points, and faces of an indexed face set.
Parameters
----------
ifs : ddg.indexedfaceset.IndexedFaceSet
coordinate_attribute : str (default="co")
Returns
-------
Mesh
Examples
--------
>>> import ddg
>>> cube = ddg.indexedfaceset.cube()
>>> points, faces, _ = ddg.arrays.from_indexed_face_set(cube)
>>> points
array([[-1, -1, -1],
[ 1, -1, -1],
[ 1, -1, 1],
[-1, -1, 1],
[-1, 1, -1],
[ 1, 1, -1],
[ 1, 1, 1],
[-1, 1, 1]])
>>> faces
array([[4, 5, 1, 0],
[6, 2, 1, 5],
[7, 3, 2, 6],
[4, 0, 3, 7],
[1, 2, 3, 0],
[7, 6, 5, 4]])
"""
# TODO: Decide if this function is even necessary.
# It shouldn't be necessary. After all, it duplicates the most basic
# functionality that one would expect from an indexed face set.
return Mesh(ifs.vertex_attributes[coordinate_key], ifs.face_list())
def _embed_points(x, /):
"""Embed points in R^d into R^3.
The embedding is::
(x,) -> (x, 0, 0) if d = 1
(x, y) -> (x, y, 0) if d = 2
(x, y, z) -> (x, y, z) if d = 3
Parameters
----------
x : numpy.ndarray of shape (d,) or (num_points, d)
Returns
-------
numpy.ndarray
"""
match x:
case np.ndarray(shape=(d,)) if d <= 3:
return np.append(x, np.zeros(3 - d))
case np.ndarray(shape=(num_points, d)) if d <= 3:
return np.hstack((x, np.zeros((num_points, 3 - d))))
case _:
raise Exception(em.should_never_happen)
def _embed(convertible, /):
match convertible:
case Points(points):
return Points(_embed_points(points))
case Curve(points, periodicity):
return Curve(_embed_points(points), periodicity)
case CurveList() as curves:
return CurveList(
[Curve(_embed_points(c.points), c.periodicity) for c in curves]
)
case Mesh(points, faces, non_manifold_edges):
return Mesh(_embed_points(points), faces, non_manifold_edges)
case _:
raise ValueError(em.should_never_happen)
# NOTE: Make sure to keep the valid input types in the docstring and the error
# message synchronised.
#
# NOTE: This has been refactored to hard-code not just the conversion
# functions, but also the sampling. The result is a much simpler signature and
# docstring.
# Previously, it was rather easy to get the sampling wrong - resulting in
# exceptions. In that case it was necessary to read another function's
# docstring just to understand how the sampling process works.
# On top of that, there are no canonical discretizations for smooth objects.
# However, most of the time it only matters that the resulting curves and
# meshes are sufficiently fine and sufficiently large in the unbounded case. If
# necessary, bounding boxes can be applied by the respective backends. If one
# needs different discretizations, one can call the other functions in this
# module, which offer more friendly function signatures and docstrings.
# If it really will be necessary to extend the signature to allow for sampling,
# it's probably better to go with
# def convert(convertible, /, *, quadric_sampling=default_quadric_sampling, etc):
# instead of
# def convert(convertible, /, *sampling)
# like before.
#
# TODO: We need to find reasonable default sampling values, which may be
# complicated. For example, for quadrics the default values might depend on the
# signature of the quadric.
#
# NOTE: This function doesn't implicitly embed R^d, where d <= 3, into R^3,
# because there may be visualization backends that can handle 2D objects out of
# the box. The new Blender conversion mechanism does embed R^d though.
[docs]def convert(
convertible,
/,
curve_sampling=_default_curve_sampling,
surface_sampling=_default_surface_sampling,
subspace_size=_default_subspace_size,
):
"""Convert pyddg objects to points, curves and meshes encoded by `array_like`.
Parameters
----------
convertible
This must be an instance of::
- `ddg.geometry.Subspace` of dimension <= 2
- `ddg.geometry.Quadric` of dimension <= 2
- `ddg.geometry.spheres.QuadricSphere` of dimension <= 2
- `ddg.indexedfaceset.IndexedFaceSet` with coordinate attribute "co"
- `ddg.halfedge.Surface` with coordinate attribute "co"
- `ddg.nets.SmoothCurve`
- `ddg.nets.SmoothNet`
- `ddg.nets.DiscreteCurve`
- `ddg.nets.DiscreteNet`
- `ddg.nets.PointNet`
- `ddg.nets.EmptyNet`
- `ddg.nets.NetCollection` of smooth or discrete nets of the same dimension
- `ddg.arrays.Points`
- `ddg.arrays.Curve`
- `ddg.arrays.CurveList`
- `ddg.arrays.Mesh`
curve_sampling : tuple of an int and a float
(default=_default_curve_sampling)
Determines the sample number and stepsize for curves
surface_sampling : tuple of an int and a float
(default=_default_surface_sampling)
Determines the sample number and stepsize for surfaces
subspace_size : float
(default=_default_subspace_size)
Determines the size of subspaces
Returns
-------
Points, Curve, CurveList or Mesh
Examples
--------
>>> import ddg
>>> line = ddg.geometry.subspace_from_affine_points((0, 0, 0), (1, 0, 0))
>>> (a, b), periodicity = ddg.arrays.convert(line)
>>> a
array([-50., 0., 0.])
>>> b
array([50., 0., 0.])
>>> periodicity
False
>>> plane = ddg.geometry.subspace_from_affine_points(
... (0, 0, 0), (1, 0, 0), (0, 1, 0)
... )
>>> points, faces, _ = ddg.arrays.convert(plane)
>>> points
array([[-50., -50., 0.],
[ 50., -50., 0.],
[ 50., 50., 0.],
[-50., 50., 0.]])
>>> faces
array([[0, 1, 2],
[0, 2, 3]])
"""
# Controls the number of samples and stepsizes used
# should this act differently on 2d objects?
N_1d = curve_sampling[0]
s_1d = curve_sampling[1]
N_2d = surface_sampling[0]
s_2d = surface_sampling[1]
L = subspace_size
match convertible:
case ddg.geometry.Subspace(dimension=0):
return from_0d_subspace(convertible)
case ddg.geometry.Subspace(dimension=1):
return from_1d_subspace(convertible, L)
case ddg.geometry.Subspace(dimension=2):
return from_2d_subspace(convertible, L, L)
case ddg.geometry.Quadric(dimension=-1):
return from_quadric(convertible)
case ddg.geometry.Quadric(dimension=0):
return from_quadric(convertible)
case ddg.geometry.Quadric(dimension=1):
return from_quadric(convertible, (N_1d, s_1d))
case ddg.geometry.Quadric(dimension=2):
return from_quadric(convertible, (N_2d, s_2d), (N_2d, s_2d))
case ddg.geometry.spheres.QuadricSphere() as sphere:
if sphere.dimension == 1:
# Scaling by sphere.radius and then rounding ensures
# that the distance between two samples is almost
# independent of sphere.radius, so that large and small
# circles will look as smooth.
num_samples = round(N_1d * sphere.radius)
else:
# For dimension 0 the number of samples doesn't matter
# and for dimension 2 we can rely on smooth shading.
num_samples = N_2d
return from_quadric_sphere(convertible, num_samples)
case ddg.indexedfaceset.IndexedFaceSet():
return from_indexed_face_set(convertible, "co")
case ddg.halfedge.Surface():
return from_half_edge_surface(convertible, "co")
case (
ddg.nets.SmoothNet()
| ddg.nets.DiscreteNet()
| ddg.nets.EmptyNet()
| ddg.nets.PointNet()
):
return _from_net(
convertible,
curve_sampling=curve_sampling,
surface_sampling=surface_sampling,
)
case ddg.nets.NetCollection():
return from_net_collection(convertible)
case Curve() | CurveList() | Mesh() | Points():
return convertible
case _:
expected = em.expected(
(
ddg.geometry.Subspace,
ddg.geometry.Quadric,
ddg.geometry.spheres.QuadricSphere,
ddg.indexedfaceset.IndexedFaceSet,
ddg.halfedge.Surface,
ddg.nets.SmoothCurve,
ddg.nets.SmoothNet,
ddg.nets.DiscreteCurve,
ddg.nets.DiscreteNet,
ddg.nets.PointNet,
ddg.nets.EmptyNet,
ddg.nets.NetCollection,
Points,
Curve,
CurveList,
Mesh,
)
)
message = em.build_error_message(
convertible, "convertible", None, [expected]
)
raise TypeError(message)
def _parse_point(point, point_name, /):
"""Parse an array representing a point in R^3.
Parameters
----------
point : array_like of shape (3,)
A point in R^3.
point_name : str
The name of point in the error message.
Returns
-------
point or (error_message : ddg._error_messages.str_)
"""
message = partial(em.build_error_message, point, point_name, em.point_definition)
point_as_array = np.asanyarray(point)
dtype = em.expected_integer_or_floating_dtype(point_as_array, point_name)
shape = em.expected_shape_at_most_3(point_as_array, point_name)
if isinstance(dtype, em.str_) or isinstance(shape, em.str_):
return message([dtype, shape])
else:
return point_as_array
def _parse_points(points, points_name, /):
"""Parse an array representing points in R^3.
Parameters
----------
points : array_like of shape (n, 3)
Points in R^3. The number n may be 0.
poinst_name : str
The name of points in the error message.
Returns
-------
points or (error_message : ddg._error_messages.str_)
"""
message = partial(em.build_error_message, points, points_name, None)
points_as_array = np.atleast_2d(np.asanyarray(points))
dtype = em.expected_integer_or_floating_dtype(points_as_array, points_name)
shape = em.expected_shape_any_at_most_3(points_as_array, points_name)
if isinstance(dtype, em.str_) or isinstance(shape, em.str_):
return message([dtype, shape])
else:
return points_as_array
def _parse_periodicity(periodicity, periodicity_name, /):
"""Parse the periodicity of a curve.
Parameters
----------
periodicity : bool
periodicity_name : str
The name of periodicity in the error message.
Returns
-------
periodicity or (error_message : ddg._error_messages.str_)
"""
return (
periodicity
if isinstance(periodicity, bool)
else em.build_error_message(
periodicity, periodicity_name, None, [em.expected(bool)]
)
)
def _parse_edges(edges, edges_name, /):
"""Parse edges.
Parameters
----------
edges : array_like of shape (n, 2) and integer dtype
Edges.
edges_name : str
The name of edges in the error message.
Returns
-------
edges or (error_message : ddg._error_messages.str_)
"""
edges_as_array = np.asanyarray(edges)
message = partial(em.build_error_message, edges_as_array, edges_name, None)
match edges_as_array:
case np.ndarray(shape=(_, 2)):
if np.issubdtype(edges_as_array.dtype, np.integer):
return edges_as_array
else:
return message([em.expected_integer_or_floating_dtype])
case np.ndarray():
return message(
[em.expected_shape_any_2, em.expected_integer_or_floating_dtype]
)
case _:
raise ValueError(em.should_never_happen)
def _is_face(sequence):
return all(isinstance(x, (int, np.integer)) for x in sequence)
def _parse_faces(faces, faces_name, /):
"""Parse faces.
Parameters
----------
faces : array_like of integer dtype or sequence of integer sequences
faces_name : str
The name of faces in the error message.
Returns
-------
faces or (error_message : ddg._error_messages.str_)
"""
# TODO: sequences with a lower-case "s" must not necessarily inherit from
# collections.abc.Sequence, see
# https://docs.python.org/3/glossary.html#term-sequence
# but we're playing somewhat loose with the types anyway...
# revisit later.
faces_as_array = np.asanyarray(faces)
if np.issubdtype(faces_as_array.dtype, np.integer):
return faces_as_array
elif isinstance(faces, Sequence) and all(_is_face(f) for f in faces):
return faces
else:
return em.build_error_message(faces, faces_name, em.faces_definition, [])
[docs]@dataclasses.dataclass
class Curve:
"""A discrete curve in R^d for d <= 3.
The curve consists of the line segments between its points.
This class should not be confused with `ddg.nets.DiscreteCurve`!
Parameters
----------
points : array_like of shape (num_points, d) and integer or floating point dtype
The dimension `d` may be at most 3.
periodicity : bool
If `True`, there is an additional line segment between the last
and the first point.
If `False`, there is no such line segment.
Attributes
----------
points : numpy.ndarray of shape (num_points, d) and integer or floating point dtype
This is an actual ndarray and not just an array_like.
periodicity : bool
"""
points: np.ndarray
periodicity: bool
def __init__(self, points, periodicity, /):
points_ = _parse_points(points, "points")
periodicity_ = _parse_periodicity(periodicity, "periodicity")
message = em.concatenate_messages(points_, periodicity_)
if isinstance(points_, em.str_) or isinstance(periodicity_, em.str_):
raise ValueError(message)
else:
self.points = points_
self.periodicity = periodicity_
def __iter__(self):
yield self.points
yield self.periodicity
[docs]class CurveList(list):
"""A list of `ddg.arrays.Curve`.
Parameters
----------
curves : sequence of ddg.array.Curve
"""
def __init__(self, curves, /):
if all(isinstance(c, Curve) for c in curves):
self.extend(curves)
else:
bad = [
(i, str(c)) for i, c in enumerate(curves) if not isinstance(c, Curve)
]
message = em.build_error_message_sequence(
bad, "curves", f"Expected instances of {Curve}"
)
raise TypeError(message)
[docs]@dataclasses.dataclass
class Mesh:
"""A mesh in R^d for d <= 3.
This is the standard representation of a polygon mesh. See
`https://en.wikipedia.org/wiki/Polygon_mesh#Face-vertex_meshes`
for more details. This class should not be confused with
`ddg.indexedfaceset.IndexedFaceSet`!
Parameters
----------
points : array_like of shape (num_points, d) and integer or floating point dtype
The dimension `d` may be at most 3.
faces : array_like or sequence
If this is an array_like, it must have shape `(num_faces, n)` and
integer dtype. In particular, every face must have n vertices. Every
row represents a face.
If this is a sequence, then its elements must themselves be sequences
of integers, which represent faces. They don't need to have the same
number of vertices.
non_manifold_edges : array_like of shape (num_edges, 2) and integer dtype
Additional edges that are distinct from the boundary edges of the faces.
Each row `(i, j)` represents an edge from the i-th to the j-th vertex.
Attributes
----------
points : numpy.ndarray of shape (num_points, d) and integer or floating point dtype
This is an actual ndarray and not just an array_like.
faces : numpy.ndarray or sequence
If this is a numpy.ndarray, then it has shape `(num_faces, n)` and
integer dtype.
If this is a sequence, then it is a sequence of integer sequences.
The numpy.ndarray representation is preferred - `faces` is always a
numpy.ndarray if every face has the same number of vertices.
non_manifold_edges : array_like of shape (num_edges, 2) and integer dtype
This is an actual ndarray and not just an array_like.
"""
points: np.ndarray
faces: np.ndarray | Sequence[Sequence[int]]
non_manifold_edges: np.ndarray = np.zeros((0, 2), int)
def __init__(self, points, faces, non_manifold_edges=np.zeros((0, 2), int), /):
maybe_points = _parse_points(points, "points")
maybe_faces = _parse_faces(faces, "faces")
maybe_edges = _parse_edges(non_manifold_edges, "non_manifold_edges")
if (
isinstance(maybe_points, em.str_)
or isinstance(maybe_faces, em.str_)
or isinstance(maybe_edges, em.str_)
):
message = em.concatenate_messages(maybe_points, maybe_faces, maybe_edges)
raise ValueError(message)
else:
self.points = maybe_points
self.faces = maybe_faces
self.non_manifold_edges = maybe_edges
def __iter__(self):
yield self.points
yield self.faces
yield self.non_manifold_edges
[docs]@dataclasses.dataclass
class Points:
"""Points in R^d for d <= 3.
Parameters
----------
points : array_like of shape (num_points, d) and integer or floating point dtype
The dimension `d` may be at most 3.
Attributes
----------
points : numpy.ndarray of shape (num_points, d) and integer or floating point dtype
This is an actual ndarray and not just an array_like.
"""
points: np.ndarray
def __init__(self, points, /):
message = partial(
em.build_error_message, points, "points", em.points_definition
)
points_as_array = np.atleast_2d(np.asanyarray(points))
dtype = em.expected_integer_or_floating_dtype(points_as_array, "points")
shape = em.expected_shape_any_at_most_3(points_as_array, "points")
if isinstance(dtype, em.str_) or isinstance(shape, em.str_):
raise ValueError(message([dtype, shape]))
else:
self.points = points_as_array
def _to_points(mesh_or_curve):
"""Converts Mesh or Curve to Points
Parameters
----------
arrays :
This must be an instance of::
- `ddg.arrays.Mesh`
- `ddg.arrays.Curve`
- `ddg.arrays.CurveList`
- `ddg.arrays.Points`
Returns
-------
ddg.arrays.Points
"""
match mesh_or_curve:
case Points():
return mesh_or_curve
case Mesh():
return Points(mesh_or_curve.points)
case Curve():
return Points(mesh_or_curve.points)
case CurveList():
return Points(np.concatenate([c.points for c in mesh_or_curve]))
case _:
expected = em.expected(
(
Points,
Curve,
CurveList,
Mesh,
)
)
message = em.build_error_message(
mesh_or_curve, "mesh_or_curve", None, [expected]
)
raise TypeError(message)
def _points_to_mesh(points, box, mesh, /):
"""Parse a point cloud.
A point mesh is a mesh representing each point.
It is put at the position of each point.
This is necessary to visualise invididual points in Blender.
Parameters
----------
points : array_like of shape (3,) or shape (_, 3)
box : ddg.arrays._BoundingBox
Only points within the box are returned. The point meshes may
stick out of the box - only their centers are guaranteed to be
within the box.
mesh : Mesh
The mesh to put at each point.
Returns
-------
(point_mesh : ddg.arrays._Mesh) or ddg._error_messages.str_
point_mesh is mesh translated by point.
"""
if not isinstance(mesh, Mesh):
raise ValueError("Invalid point mesh for PointCloud Visual")
if not isinstance(points, Points):
raise ValueError(f"Expected a Points not a {type(points)}")
else:
if box.is_R3:
bounded_points = points.points
else:
bounded_points = _bounding_box_points(points.points, box)
meshes = [
Mesh(mesh.points + p, mesh.faces, mesh.non_manifold_edges)
for p in bounded_points
]
point_cloud = mesh_union(meshes)
if isinstance(point_cloud, Mesh):
return point_cloud
else:
raise Exception(em.should_never_happen)
[docs]def mesh_union(meshes, /):
"""Union of meshes.
Parameters
----------
meshes : sequence of ddg.arrays.Mesh
Returns
-------
ddg.arrays.Mesh
Examples
--------
>>> import ddg
>>> triangle_1 = ddg.arrays.Mesh(
... np.array([(0, 0, 0), (1, 0, 0), (1, 1, 0)]), [(0, 1, 2)]
... )
>>> triangle_2 = ddg.arrays.Mesh(
... np.array([(0, 0, 0), (0, 1, 0), (0, 0, 1)]), [(0, 1, 2)]
... )
>>> points, faces, non_manifold_edges = ddg.arrays.mesh_union(
... (triangle_1, triangle_2)
... )
>>> points
array([[0, 0, 0],
[1, 0, 0],
[1, 1, 0],
[0, 0, 0],
[0, 1, 0],
[0, 0, 1]])
>>> faces
array([[0, 1, 2],
[3, 4, 5]])
>>> non_manifold_edges
array([], shape=(0, 2), dtype=int64)
"""
if not all(isinstance(m, Mesh) for m in meshes):
bad = [(i, m) for i, m in enumerate(meshes) if not isinstance(m, tuple)]
message = em.build_error_message_sequence(
bad, "meshes", f"Expected instances of {Mesh}"
)
raise TypeError(message)
elif meshes:
points = np.vstack([mesh.points for mesh in meshes])
# We need to fix the indices of every edge and face not part of the
# first mesh because the corresponding points are no longer indexed
# from 0.
num_points = np.cumsum([mesh.points.shape[0] for mesh in meshes])
shifts = np.insert(num_points[:-1], 0, 0)
shifted_edges = [
mesh.non_manifold_edges + n for mesh, n in zip(meshes, shifts, strict=True)
]
edges = np.vstack(shifted_edges)
shifted_faces = [
mesh.faces + n
if isinstance(mesh.faces, np.ndarray)
else tuple(v + n for v in mesh.faces)
for mesh, n in zip(meshes, shifts, strict=True)
]
def non_empty(faces):
return (isinstance(faces, np.ndarray) and faces.shape[0] > 0) or faces != ()
non_empty_shifted_faces = list(filter(non_empty, shifted_faces))
try:
# The good case where every face has the same number of
# vertices/edges, so we can just vstack everything.
faces = np.vstack(non_empty_shifted_faces)
except ValueError:
# The bad case: if there are two faces with different
# numbers of vertices (respectively edges), vstack fails
# with
#
# ValueError: all the input array dimensions for the concatenation
# axis must match exactly, [...]
#
# We have to fall back to the tuple representation.
tuple_faces = [tuple(tuple(f) for f in fs) for fs in shifted_faces]
faces = sum(tuple_faces, ())
return Mesh(points, faces, edges)
else:
# If meshes is empty, return an empty mesh in R^3. Because vstack
# cannot stack 0 arrays, we need to handle this separately.
return Mesh(np.zeros((0, 3), float), (), np.zeros((0, 2), int))
[docs]def edges(mesh, /):
non_unique_face_edges = np.vstack(
[np.column_stack([f, np.roll(f, -1)]) for f in mesh.faces]
)
unique_face_edges = np.unique(non_unique_face_edges, axis=0)
return np.vstack([unique_face_edges, mesh.non_manifold_edges])
@dataclasses.dataclass
class _BoundingBox:
min_: np.ndarray
max_: np.ndarray
@property
def center(self):
return np.nan_to_num(self.min_ + (self.max_ - self.min_) / 2)
@property
def distances(self):
return self.max_ - self.center
@property
def is_R3(self):
return np.all(self.min_ == -np.inf) and np.all(self.max_ == np.inf)
def _bounding_box_points(points, box, /):
"""Filter the the subset of points that lies in box.
Parameters
----------
points : numpy.ndarray of shape (_, 3) and integer or floating point dtype
box : ddg.arrays._BoundingBox
Returns
-------
numpy.ndarray of shape (_, 3) and integer or floating point dtype
"""
# There are more concise and probably faster vectorised
# implementations, but this is easy to understand.
x, y, z = points.T
x_min, y_min, z_min = box.min_
x_max, y_max, z_max = box.max_
x_good = (x_min <= x) & (x <= x_max)
y_good = (y_min <= y) & (y <= y_max)
z_good = (z_min <= z) & (y <= z_max)
good_mask = x_good & y_good & z_good
return points[good_mask]
def _bounding_box_segment_1d(a, b, min_, max_, /):
"""Solve min_ <= a + t * (b - a) <= max_ for t.
Parameters
----------
a, b, min_, max_ : float
Returns
-------
tuple[float, float] or None
If there is a solution, return the solution interval as pair of floats.
Otherwise return None.
"""
if np.isclose(a, b):
if min_ <= a <= max_:
return -np.inf, np.inf
else:
return None
else:
d = b - a
if d >= 0:
t_min = (min_ - a) / d
t_max = (max_ - a) / d
else:
t_min = (max_ - a) / d
t_max = (min_ - a) / d
return t_min, t_max
def _bounding_box_segment(a, b, box, /):
"""Compute the intersection of the segment ab with box.
Parameters
----------
a : numpy.ndarray of shape (3,)
b : numpy.ndarray of shape (3,)
box : ddg.arrays._BoundingBox
Returns
-------
tuple[numpy.ndarray, numpy.ndarray] or None
If the segment intersects the box in a subsegment, return the
subsegment p, q such that a, p, q, b lie in that order on the line
defined by ab.
Otherwise return None.
"""
# Consider the line defined by ab. We wish to intersect it with the
# bounding box.
# The line is parametrised by t -> a + t * (b - a), so this amounts to
# solving
#
# box.min_[i] <= a[i] + t_i * (b[i] - a[i]) <= box.max_[i]
#
# where i = 0, 1, 2.
# It is possible that some of these inequalities have no solution. This
# happens if and only if if b[i] == a[i] not in [box.min_[i], box.max_[i]]
# for some i.
intervals = [
_bounding_box_segment_1d(a[i], b[i], box.min_[i], box.max_[i]) for i in range(3)
]
if None in intervals:
return None
else:
# The above inequalities have solutions for all i.
# Their intersection of the three solution intervals is [t_min, t_max].
t_min = np.max([t for t, _ in intervals])
t_max = np.min([t for _, t in intervals])
if t_min > t_max:
# The common solution interval is empty, so the line ab doesn't
# intersect the box.
return None
else:
# The line ab intersects the box. It remains to check if the
# segment (a, b) intersects the box.
# The parametrisation
# R -> line ab, t -> a + t * (b - a)
# maps 0 to a and 1 to b. The comments below visualise the order of 0,
# 1, t_min, t_max on R in the respective cases.
if 1 < t_min:
# 0 1 t_min t_max
return None
elif 0 > t_max:
# t_min t_max 0 1
return None
else:
# Covers all four remaining cases:
# 0 t_min t_max 1
# 0 t_min 1 t_max
# t_min 0 t_max 1
# t_min 0 1 t_max
d = b - a
t_min_ = max(t_min, 0)
t_max_ = min(t_max, 1)
return a + t_min_ * d, a + t_max_ * d
def _components(segments_, /):
"""Combine consecutive segments into a curve.
Parameters
----------
segments_ : sequence[tuple[numpy.ndarray, np.ndarray]]
Segments in R^3.
Returns
-------
numpy.ndarray of shape (n, 3) and sequence[tuple[numpy.ndarray, np.ndarray]]
The first value is a curve comprising of the union of the first n + 1
segments of segments_. The remaining segments are the second value.
"""
match segments_:
case [(a, b), *remaining]:
points = [a, b]
for (c, d) in remaining:
last = points[-1]
if np.allclose(last, c):
# curve continues
points.append(d)
else:
# curve ends
break
return np.array(points), remaining
case _:
raise Exception(em.should_never_happen)
def _bounding_box_Curve(curve, box, /):
"""Intersect a curve with a bounding box.
Parameters
----------
curve : ddg.arrays.Curve
box : ddg.arrays._BoundingBox
Returns
-------
ddg.arrays.DisconnectedCurve
"""
# The algorithm and implementation is not efficient, but it's relatively
# simple and it works.
# Basic idea:
# A curve comprises of line segments.
# - Compute all the subsegments intersecting the box.
# - Combine adjacent segments if their end and start are sufficiently close.
# - Otherwise adjacent segments belong to different components of the
# intersection.
points = curve.points
maybe_segments = [
_bounding_box_segment(a, b, box) for a, b in zip(points[:-1], points[1:])
]
segments = [s for s in maybe_segments if s is not None]
components = []
while segments:
c, remaining_segments = _components(segments)
components.append(c)
segments = remaining_segments
match components:
case []:
return CurveList([])
case [c]:
# One component means that every segment p_0p_1, ...,
# p_{n-2}p_{n-1} is contained in box. Since box is convex, that
# means that p_{n-1}p_0 is also contained in the box. This means
# that the curve is contained in the box. We return curve directly,
# so the periodicity is correct and we avoid rounding errors
# accrued in the computation of components.
return CurveList([Curve(c, False)])
case _:
if curve.periodicity:
# c and c_ are distinct because len(components) >= 2 in this branch.
c = curve[0]
c_ = curve[-1]
if np.allclose(c[0], c_[:-1]):
# c and c_ are actually two parts of one component of curve
# intersected with box because curve is periodic.
components_ = [c_ + c] + components[1:-1]
else:
components_ = components
return CurveList([Curve(c, False) for c in components_])
else:
# No component can be periodic.
return CurveList([Curve(c, False) for c in components])
def _bounding_box_Curve_DisconnectedCurve(curve, box, /):
"""Intersect a curve or disconnected curve with a bounding box.
Parameters
----------
curve : ddg.arrays.Curve | ddg.arrays.DisconnectedCurve
box : ddg.arrays._BoundingBox
Returns
-------
ddg.arrays.DisconnectedCurve
"""
match curve:
case Curve() as curve:
return _bounding_box_Curve(curve, box)
case CurveList() as curves:
components = [_bounding_box_Curve(curve, box) for curve in curves]
return CurveList(sum(components, []))
case _:
raise ValueError(em.should_never_happen)
[docs]def line_segment_from_points(point0, point1, /, domain=(0, 1)):
"""Get a line segment defined by two points as a curve.
Returns the line segment between the two points::
(1 - domain[0]) * point0 + domain[0] * point1,
(1 - domain[1]) * point0 + domain[1] * point1.
Parameters
----------
point0, point1 : ddg.geometry.Point, Points, or array_like
The points defining the line segment or their affine coordinates.
Must be at most 3-dimensional.
domain : tuple[float, float], optinal
Defines the end points of the line segment in relation to the first point and
the direction vector between the two points.
Returns
-------
ddg.arrays.Curve
Raises
------
ValueError
If a given point is at infinity.
Examples
--------
>>> import ddg
>>> line_segment = ddg.arrays.line_segment_from_points(
... [0, 0, 0], ddg.geometry.Point([2, 2, 2, 1]), [-0.5, 1]
... )
>>> line_segment
Curve(points=array([[-1., -1., -1.],
[ 2., 2., 2.]]), periodicity=False)
"""
if isinstance(point0, ddg.geometry.Point):
if point0.at_infinity():
raise ValueError(
"A point at infinity was given. Consider using the function"
"line_segment_from_point_and_direction."
)
point0_coords = point0.affine_point
elif isinstance(point0, Points) and point0.points.shape[0] == 1:
point0_coords = point0.points[0]
else:
point0_coords = _parse_point(point0, "point0")
if isinstance(point0_coords, em.str_):
raise ValueError(point0_coords)
if isinstance(point1, ddg.geometry.Point):
if point1.at_infinity():
raise ValueError(
"A point at infinity was given. Consider using the function"
"line_segment_from_point_and_direction."
)
point1_coords = point1.affine_point
elif isinstance(point1, Points) and point1.points.shape[0] == 1:
point1_coords = point1.points[0]
else:
point1_coords = _parse_point(point1, "point1")
if isinstance(point0_coords, em.str_):
raise ValueError(point0_coords)
direction = point1_coords - point0_coords
curve = Curve(
(
point0_coords + domain[0] * direction,
point0_coords + domain[1] * direction,
),
False,
)
return curve
[docs]def line_segment_from_point_and_direction(
point, direction, domain=(0, 1), normalized=False
):
"""Get a line segment defined by a point and a direction vector as a curve.
Returns the line segment between the two points::
point + domain[0] * direction,
point + domain[1] * direction
where `direction` may have been normalized before.
Parameters
----------
point : ddg.geometry.Point, Points or array_like
Location point of the line segment. Must be at most 3-dimensional.
direction : ddg.geometry.Point or array_like
Direction vector of the line segment as a `Point` at infinity or as a
coordinate vector. The dimension must match the one of `point`.
domain : tuple[float, float], optional
Defines the end points of the line segment in relation to the location point
and the direction vector.
normalized : bool, optional
Whether the direction vector should be normalized. If this is set to ``True``,
the length of the line segment coincides with the length of the domain.
Returns
-------
ddg.arrays.Curve
Raises
------
ValueError
If the location point lies at infinity.
If the direction vector is an object of type `ddg.geometry.Point` which is
not at infinity.
Examples
--------
>>> import ddg
>>> line_segment = ddg.arrays.line_segment_from_point_and_direction(
... [1, 1, 1], ddg.geometry.Point([1, 1, 1, 0]), [-1, 1]
... )
>>> line_segment
Curve(points=array([[0., 0., 0.],
[2., 2., 2.]]), periodicity=False)
>>> line_segment = ddg.arrays.line_segment_from_point_and_direction(
... [1, 1, 1], ddg.geometry.Point([10, 10, 10, 0]), [-1, 1], True
... )
>>> line_segment
Curve(points=array([[0.42264973, 0.42264973, 0.42264973],
[1.57735027, 1.57735027, 1.57735027]]), periodicity=False)
"""
if isinstance(point, ddg.geometry.Point):
if point.at_infinity():
raise ValueError("point must not be at infinity.")
point_coords = point.affine_point
elif isinstance(point, Points) and point.points.shape[0] == 1:
point_coords = point.points[0]
else:
point_coords = _parse_point(point, "point")
if isinstance(point_coords, em.str_):
raise ValueError(point_coords)
if isinstance(direction, ddg.geometry.Point):
if not direction.at_infinity():
raise ValueError("direction must be a point at infinity")
direction_vector = direction.point[:-1]
else:
direction_vector = _parse_point(direction, "direction")
if isinstance(direction_vector, em.str_):
raise ValueError(direction_vector)
if normalized:
direction_vector /= np.linalg.norm(direction_vector)
curve = Curve(
(
point_coords + domain[0] * direction_vector,
point_coords + domain[1] * direction_vector,
),
False,
)
return curve