from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
from math import cos
from math import pi
from copy import deepcopy
from compas.geometry import centroid_points
from compas.geometry import offset_polygon
from compas.utilities import iterable_like
from compas.utilities import pairwise
from .core import BaseMesh
__all__ = [
'mesh_subdivide',
'mesh_subdivide_tri',
'mesh_subdivide_corner',
'mesh_subdivide_quad',
'mesh_subdivide_catmullclark',
'mesh_subdivide_doosabin',
'mesh_subdivide_frames',
'trimesh_subdivide_loop',
]
def mesh_fast_copy(other):
subd = SubdMesh()
subd.vertex = deepcopy(other.vertex)
subd.face = deepcopy(other.face)
# subd.edgedata = deepcopy(other.edgedata)
subd.facedata = deepcopy(other.facedata)
subd.halfedge = deepcopy(other.halfedge)
subd._max_vertex = other._max_vertex
subd._max_face = other._max_face
return subd
class SubdMesh(BaseMesh):
from .core import mesh_split_edge
_add_vertex = BaseMesh.add_vertex
_add_face = BaseMesh.add_face
_insert_vertex = BaseMesh.insert_vertex
split_edge = mesh_split_edge
def add_vertex(self, x, y, z):
key = self._max_vertex = self._max_vertex + 1
if key not in self.vertex:
self.vertex[key] = {}
self.halfedge[key] = {}
self.vertex[key] = dict(x=x, y=y, z=z)
return key
def add_face(self, vertices):
fkey = self._max_face = self._max_face + 1
self.face[fkey] = vertices
self.facedata[fkey] = {}
for i in range(-1, len(vertices) - 1):
u = vertices[i]
v = vertices[i + 1]
self.halfedge[u][v] = fkey
if u not in self.halfedge[v]:
self.halfedge[v][u] = None
return fkey
def insert_vertex(self, fkey):
x, y, z = self.face_center(fkey)
w = self.add_vertex(x=x, y=y, z=z)
for u, v in self.face_halfedges(fkey):
self.add_face([u, v, w])
del self.face[fkey]
return w
# distinguish between subd of meshes with and without boundary
# closed vs. open
# pay attention to extraordinary points
# and to special rules on boundaries
# interpolation vs. approxmation?!
# add numerical versions to compas.datastructures.mesh.(algorithms.)numerical
# investigate meaning and definition of limit surface
# any subd algorithm should return a new subd mesh, leaving the control mesh intact
[docs]def mesh_subdivide(mesh, scheme='catmullclark', **options):
"""Subdivide the input mesh.
Parameters
----------
mesh : Mesh
A mesh object.
scheme : {'tri', 'corner', 'catmullclark', 'doosabin'}, optional
The scheme according to which the mesh should be subdivided.
Default is ``'catmullclark'``.
options : dict
Optional additional keyword arguments.
Returns
-------
Mesh
The subdivided mesh.
Raises
------
NotImplementedError
If the scheme is not supported.
"""
if scheme == 'tri':
return mesh_subdivide_tri(mesh, **options)
if scheme == 'quad':
return mesh_subdivide_quad(mesh, **options)
if scheme == 'corner':
return mesh_subdivide_corner(mesh, **options)
if scheme == 'catmullclark':
return mesh_subdivide_catmullclark(mesh, **options)
if scheme == 'doosabin':
return mesh_subdivide_doosabin(mesh, **options)
raise NotImplementedError
[docs]def mesh_subdivide_tri(mesh, k=1):
"""Subdivide a mesh using simple insertion of vertices.
Parameters
----------
mesh : Mesh
The mesh object that will be subdivided.
k : int
Optional. The number of levels of subdivision. Default is ``1``.
Returns
-------
Mesh
A new subdivided mesh.
Examples
--------
>>> box = Box.from_corner_corner_height([0.0, 0.0, 0.0], [1.0, 1.0, 0.0], 1.0)
>>> mesh = Mesh.from_shape(box)
>>> k = 2
>>> subd = mesh_subdivide_tri(mesh, k=k)
>>> mesh is subd
False
>>> type(mesh) is type(subd)
True
>>> k1 = sum(len(mesh.face_vertices(fkey)) for fkey in mesh.faces())
>>> subd.number_of_faces() == (k1 if k == 1 else k1 * 3 ** (k - 1))
True
"""
cls = type(mesh)
subd = mesh_fast_copy(mesh)
for _ in range(k):
for fkey in list(subd.faces()):
subd.insert_vertex(fkey)
return cls.from_data(subd.data)
[docs]def mesh_subdivide_quad(mesh, k=1):
"""Subdivide a mesh such that all faces are quads.
Parameters
----------
mesh : Mesh
The mesh object that will be subdivided.
k : int
Optional. The number of levels of subdivision. Default is ``1``.
Returns
-------
Mesh
A new subdivided mesh.
Examples
--------
>>> box = Box.from_corner_corner_height([0.0, 0.0, 0.0], [1.0, 1.0, 0.0], 1.0)
>>> mesh = Mesh.from_shape(box)
>>> k = 2
>>> subd = mesh_subdivide_quad(mesh, k=k)
>>> mesh is subd
False
>>> type(mesh) is type(subd)
True
>>> subd.number_of_faces() == mesh.number_of_faces() * 4 ** k
True
"""
cls = type(mesh)
subd = mesh_fast_copy(mesh)
for face in subd.faces():
subd.facedata[face]['path'] = [face]
for _ in range(k):
faces = {face: subd.face_vertices(face)[:] for face in subd.faces()}
face_centroid = {face: subd.face_centroid(face) for face in subd.faces()}
for u, v in list(subd.edges()):
subd.split_edge(u, v, allow_boundary=True)
for face, vertices in faces.items():
descendant = {i: j for i, j in subd.face_halfedges(face)}
ancestor = {j: i for i, j in subd.face_halfedges(face)}
x, y, z = face_centroid[face]
c = subd.add_vertex(x=x, y=y, z=z)
for i, vertex in enumerate(vertices):
a = ancestor[vertex]
d = descendant[vertex]
newface = subd.add_face([a, vertex, d, c])
subd.facedata[newface]['path'] = subd.facedata[face]['path'] + [i]
del subd.face[face]
del subd.facedata[face]
subd2 = cls.from_data(subd.data)
return subd2
[docs]def mesh_subdivide_corner(mesh, k=1):
"""Subdivide a mesh by cutting corners.
Parameters
----------
mesh : Mesh
The mesh object that will be subdivided.
k : int
Optional. The number of levels of subdivision. Default is ``1``.
Returns
-------
Mesh
A new subdivided mesh.
Notes
-----
This is essentially the same as Loop subdivision, but applied to general
meshes.
"""
cls = type(mesh)
for _ in range(k):
subd = mesh_fast_copy(mesh)
# split every edge
for u, v in list(subd.edges()):
subd.split_edge(u, v, allow_boundary=True)
# create 4 new faces for every old face
for fkey in mesh.faces():
descendant = {i: j for i, j in subd.face_halfedges(fkey)}
ancestor = {j: i for i, j in subd.face_halfedges(fkey)}
center = []
for key in mesh.face_vertices(fkey):
a = ancestor[key]
d = descendant[key]
subd.add_face([a, key, d])
center.append(a)
subd.add_face(center)
del subd.face[fkey]
mesh = subd
subd2 = cls.from_data(mesh.data)
return subd2
[docs]def mesh_subdivide_catmullclark(mesh, k=1, fixed=None):
"""Subdivide a mesh using the Catmull-Clark algorithm.
Parameters
----------
mesh : Mesh
The mesh object that will be subdivided.
k : int
Optional. The number of levels of subdivision. Default is ``1``.
fixed : list
Optional. A list of fixed vertices. Default is ``None``.
Returns
-------
Mesh
A new subdivided mesh.
Notes
-----
Note that *Catmull-Clark* subdivision is like *Quad* subdivision, but with
smoothing after every level of further subdivision. Smoothing is done
according to the scheme prescribed by the Catmull-Clark algorithm.
Examples
--------
>>> box = Box.from_corner_corner_height([0.0, 0.0, 0.0], [1.0, 1.0, 0.0], 1.0)
>>> mesh = Mesh.from_shape(box)
>>> k = 2
>>> subd = mesh_subdivide_catmullclark(mesh, k=k)
>>> mesh is subd
False
>>> type(mesh) is type(subd)
True
>>> subd.number_of_faces() == mesh.number_of_faces() * 4 ** k
True
"""
cls = type(mesh)
if not fixed:
fixed = []
fixed = set(fixed)
for _ in range(k):
subd = mesh_fast_copy(mesh)
# keep track of original connectivity and vertex locations
bkeys = set(subd.vertices_on_boundary())
bkey_edgepoints = {key: [] for key in bkeys}
# apply quad meshivision scheme
# keep track of the created edge points that are not on the boundary
# keep track track of the new edge points on the boundary
# and their relation to the previous boundary points
# quad subdivision
# ======================================================================
edgepoints = []
for u, v in mesh.edges():
w = subd.split_edge(u, v, allow_boundary=True)
# document why this is necessary
# everything else in this loop is just quad subdivision
if u in bkeys and v in bkeys:
bkey_edgepoints[u].append(w)
bkey_edgepoints[v].append(w)
continue
edgepoints.append(w)
fkey_xyz = {fkey: mesh.face_centroid(fkey) for fkey in mesh.faces()}
for fkey in mesh.faces():
descendant = {i: j for i, j in subd.face_halfedges(fkey)}
ancestor = {j: i for i, j in subd.face_halfedges(fkey)}
x, y, z = fkey_xyz[fkey]
c = subd.add_vertex(x=x, y=y, z=z)
for key in mesh.face_vertices(fkey):
a = ancestor[key]
d = descendant[key]
subd.add_face([a, key, d, c])
del subd.face[fkey]
# update coordinates
# ======================================================================
# these are the coordinates before updating
key_xyz = {key: subd.vertex_coordinates(key) for key in subd.vertex}
# move each edge point to the average of the neighboring centroids and
# the original end points
for w in edgepoints:
x, y, z = centroid_points([key_xyz[nbr] for nbr in subd.halfedge[w]])
subd.vertex[w]['x'] = x
subd.vertex[w]['y'] = y
subd.vertex[w]['z'] = z
# move each vertex to the weighted average of itself, the neighboring
# centroids and the neighboring mipoints
for key in mesh.vertices():
if key in fixed:
continue
if key in bkeys:
nbrs = set(bkey_edgepoints[key])
nbrs = [key_xyz[nbr] for nbr in nbrs]
e = 0.5
v = 0.5
E = [coord * e for coord in centroid_points(nbrs)]
V = [coord * v for coord in key_xyz[key]]
x, y, z = [E[_] + V[_] for _ in range(3)]
else:
fnbrs = [mesh.face_centroid(fkey) for fkey in mesh.vertex_faces(key) if fkey is not None]
nbrs = [key_xyz[nbr] for nbr in subd.halfedge[key]]
n = float(len(nbrs))
f = 1.0 / n
e = 2.0 / n
v = (n - 3.0) / n
F = centroid_points(fnbrs)
E = centroid_points(nbrs)
V = key_xyz[key]
x = f * F[0] + e * E[0] + v * V[0]
y = f * F[1] + e * E[1] + v * V[1]
z = f * F[2] + e * E[2] + v * V[2]
subd.vertex[key]['x'] = x
subd.vertex[key]['y'] = y
subd.vertex[key]['z'] = z
mesh = subd
subd2 = cls.from_data(mesh.data)
return subd2
[docs]def mesh_subdivide_doosabin(mesh, k=1, fixed=None):
"""Subdivide a mesh following the doo-sabin scheme.
Parameters
----------
mesh : Mesh
The mesh object that will be subdivided.
k : int
Optional. The number of levels of subdivision. Default is ``1``.
fixed : list
Optional. A list of fixed vertices. Default is ``None``.
Returns
-------
Mesh
A new subdivided mesh.
Examples
--------
>>> box = Box.from_corner_corner_height([0.0, 0.0, 0.0], [1.0, 1.0, 0.0], 1.0)
>>> mesh = Mesh.from_shape(box)
>>> k = 2
>>> subd = mesh_subdivide_doosabin(mesh, k=k)
>>> mesh is subd
False
>>> type(mesh) is type(subd)
True
"""
if not fixed:
fixed = []
fixed = set(fixed)
cls = type(mesh)
for _ in range(k):
old_xyz = {key: mesh.vertex_coordinates(key) for key in mesh.vertices()}
fkey_old_new = {fkey: {} for fkey in mesh.faces()}
subd = SubdMesh()
for fkey in mesh.faces():
vertices = mesh.face_vertices(fkey)
n = len(vertices)
face = []
for i in range(n):
old = vertices[i]
cx, cy, cz = 0, 0, 0
for j in range(n):
x, y, z = old_xyz[vertices[j]]
if i == j:
alpha = (n + 5.) / (4. * n)
else:
alpha = (3. + 2. * cos(2. * pi * (i - j) / n)) / (4. * n)
cx += alpha * x
cy += alpha * y
cz += alpha * z
new = subd.add_vertex(cx, cy, cz)
fkey_old_new[fkey][old] = new
face.append(new)
subd.add_face(face)
boundary = set(mesh.vertices_on_boundary())
for key in mesh.vertices():
if key in boundary:
continue
face = [fkey_old_new[fkey][key] for fkey in mesh.vertex_faces(key, ordered=True) if fkey is not None]
subd.add_face(face[::-1])
edges = set()
for u in mesh.halfedge:
for v in mesh.halfedge[u]:
if (u, v) in edges:
continue
edges.add((u, v))
edges.add((v, u))
uv_fkey = mesh.halfedge[u][v]
vu_fkey = mesh.halfedge[v][u]
if uv_fkey is None or vu_fkey is None:
continue
face = [
fkey_old_new[uv_fkey][u],
fkey_old_new[vu_fkey][u],
fkey_old_new[vu_fkey][v],
fkey_old_new[uv_fkey][v]
]
subd.add_face(face)
mesh = subd
subd2 = cls.from_data(mesh.data)
return subd2
def mesh_subdivide_frames(mesh, offset, add_windows=False):
"""Subdivide a mesh by creating offset frames and windows on its faces.
Parameters
----------
mesh : Mesh
The mesh object to be subdivided.
offset : float or dict
The offset distance to create the frames.
A single value will result in a constant offset everywhere.
A dictionary mapping facekey: offset will be processed accordingly.
add_windows : boolean
Optional. Flag to add window face. Default is ``False``.
Returns
-------
Mesh
A new subdivided mesh.
Examples
--------
>>>
"""
subd = SubdMesh()
# 0. pre-compute offset distances
if not isinstance(offset, dict):
distances = iterable_like(mesh.faces(), [offset], offset)
offset = {fkey: od for fkey, od in zip(mesh.faces(), distances)}
# 1. add vertices
newkeys = {}
for vkey, attr in mesh.vertices(True):
newkeys[vkey] = subd.add_vertex(*mesh.vertex_coordinates(vkey))
# 2. add faces
for fkey in mesh.faces():
face = [newkeys[vkey] for vkey in mesh.face_vertices(fkey)]
d = offset.get(fkey)
# 2a. add face and break if no offset is found
if d is None:
subd.add_face(face)
continue
polygon = offset_polygon(mesh.face_coordinates(fkey), d)
# 2a. add offset vertices
window = []
for xyz in polygon:
x, y, z = xyz
new_vkey = subd.add_vertex(x=x, y=y, z=z)
window.append(new_vkey)
# 2b. frame faces
face = face + face[:1]
window = window + window[:1]
for sa, sb in zip(pairwise(face), pairwise(window)):
subd.add_face([sa[0], sa[1], sb[1], sb[0]])
# 2c. window face
if add_windows:
subd.add_face(window)
return subd
def trimesh_subdivide_loop(mesh, k=1, fixed=None):
"""Subdivide a triangle mesh using the Loop algorithm.
Parameters
----------
mesh : Mesh
The mesh object that will be subdivided.
k : int
Optional. The number of levels of subdivision. Default is ``1``.
fixed : list
Optional. A list of fixed vertices. Default is ``None``.
Returns
-------
Mesh
A new subdivided mesh.
Examples
--------
Make a low poly mesh from a box shape.
Triangulate the faces.
>>> box = Box.from_corner_corner_height([0.0, 0.0, 0.0], [1.0, 1.0, 0.0], 1.0)
>>> mesh = Mesh.from_shape(box)
>>> mesh_quads_to_triangles(mesh)
Subdivide 2 times.
>>> k = 2
>>> subd = trimesh_subdivide_loop(mesh, k=k)
Compare low-poly cage with subdivision mesh.
>>> mesh is subd
False
>>> type(mesh) is type(subd)
True
>>> subd.number_of_faces() == mesh.number_of_faces() * (3 + 1) ** k
True
"""
cls = type(mesh)
if not fixed:
fixed = []
fixed = set(fixed)
subd = mesh_fast_copy(mesh)
for _ in range(k):
key_xyz = {key: subd.vertex_coordinates(key) for key in subd.vertices()}
fkey_vertices = {fkey: subd.face_vertices(fkey)[:] for fkey in subd.faces()}
uv_w = {(u, v): subd.face_vertex_ancestor(fkey, u) for fkey in subd.faces() for u, v in subd.face_halfedges(fkey)}
boundary = set(subd.vertices_on_boundary())
for key in subd.vertices():
nbrs = subd.vertex_neighbors(key)
if key in boundary:
xyz = key_xyz[key]
x = 0.75 * xyz[0]
y = 0.75 * xyz[1]
z = 0.75 * xyz[2]
for n in nbrs:
if subd.halfedge[key][n] is None or subd.halfedge[n][key] is None:
xyz = key_xyz[n]
x += 0.125 * xyz[0]
y += 0.125 * xyz[1]
z += 0.125 * xyz[2]
else:
n = len(nbrs)
if n == 3:
a = 3. / 16.
else:
a = 3. / (8 * n)
xyz = key_xyz[key]
nbrs = [key_xyz[nbr] for nbr in nbrs]
nbrs = [sum(axis) for axis in zip(*nbrs)]
x = (1. - n * a) * xyz[0] + a * nbrs[0]
y = (1. - n * a) * xyz[1] + a * nbrs[1]
z = (1. - n * a) * xyz[2] + a * nbrs[2]
subd.vertex[key]['x'] = x
subd.vertex[key]['y'] = y
subd.vertex[key]['z'] = z
edgepoints = {}
# odd vertices
for u, v in list(subd.edges()):
w = subd.split_edge(u, v, allow_boundary=True)
edgepoints[(u, v)] = w
edgepoints[(v, u)] = w
a = key_xyz[u]
b = key_xyz[v]
if (u, v) in uv_w and (v, u) in uv_w:
c = key_xyz[uv_w[(u, v)]]
d = key_xyz[uv_w[(v, u)]]
xyz = [(3.0 / 8.0) * (a[i] + b[i]) + (1.0 / 8.0) * (c[i] + d[i]) for i in range(3)]
else:
xyz = [0.5 * (a[i] + b[i]) for i in range(3)]
subd.vertex[w]['x'] = xyz[0]
subd.vertex[w]['y'] = xyz[1]
subd.vertex[w]['z'] = xyz[2]
# new faces
for fkey, vertices in fkey_vertices.items():
u, v, w = vertices
uv = edgepoints[(u, v)]
vw = edgepoints[(v, w)]
wu = edgepoints[(w, u)]
subd.add_face([wu, u, uv])
subd.add_face([uv, v, vw])
subd.add_face([vw, w, wu])
subd.add_face([uv, vw, wu])
del subd.face[fkey]
subd2 = cls.from_data(subd.data)
return subd2
# ==============================================================================
# Main
# ==============================================================================
if __name__ == "__main__":
import doctest
import compas # noqa: F401
from compas.datastructures import mesh_quads_to_triangles # noqa: F401
from compas.datastructures import Mesh # noqa: F401
from compas.geometry import Box # noqa: F401
doctest.testmod(globs=globals())
# from compas.datastructures import Mesh
# from compas.geometry import Box
# from compas.utilities import print_profile
# from compas_viewers.multimeshviewer import MultiMeshViewer
# subdivide = print_profile(mesh_subdivide_quad)
# box = Box.from_width_height_depth(10.0, 10.0, 10.0)
# mesh = Mesh.from_shape(box)
# mesh.default_face_attributes.update(path=[])
# subd = subdivide(mesh, k=3)
# print(subd.face_attribute(subd.get_any_face(), 'path'))
# # print(mesh.number_of_faces())
# viewer = MultiMeshViewer()
# viewer.meshes = [subd]
# viewer.show()