"""Module `delaunay` gives methods at hand to check for Delaunay property of
a triangulation and to transform it into a Delaunay triangulation
via edge-flips. """
import logging
import math
import os
# import argparse
from pathlib import Path
from warnings import warn
import numpy as np
import ddg
import ddg.halfedge._copy
import ddg.halfedge._io as io
import ddg.halfedge._modify
import ddg.halfedge._set
from ddg.halfedge._surface import Surface
from ddg.math import euclidean2d
logger = logging.getLogger(__name__)
[docs]def calculate_angle(e):
"""
Calculates the inner angle between an edge and its consecutive next edge
within a face.
Parameters
----------
e : edge
The edge pointing towards the inner angle that will be calculated.
Returns
-------
angle : double
The angle between `e` and `e.nex`.
Raises
------
ValueError
If the edge is a boundary edge, i.e. does not belong to a face.
If the triangle is degenerated or the triangle inequality doesn't hold.
"""
if e.face is None:
raise ValueError(
"The edge ({} {}) does not belong to a face, thus no angle "
"will be calculated.".format(e.tail, e.head)
)
if not hasattr(e, "len"):
raise ValueError(
"The edge ({} {}) does not have the attribute `len` "
"which is necessary for calculation.".format(e.tail, e.head)
)
a = e.len
b = e.nex.len
c = e.pre.len
if a * b == 0:
raise ValueError(
"The triangle ({} {} {}) is degenerated.".format(
e.pre.head, e.head, e.nex.head
)
)
if abs((a * a + b * b - c * c) / (2 * a * b)) > 1:
raise ValueError(
"The triangle inequality does not hold for the "
"triangle with edges ({}{}, {}{}, {}{}).".format(
e.pre.tail, e.pre.head, e.tail, e.head, e.nex.tail, e.nex.head
)
)
angle = abs(math.acos((a**2 + b**2 - c**2) / (2 * a * b)))
return angle
[docs]def opp_halfangle_tan(e):
r"""
Calculates the tangent of half the angle that lies opposite the given
halfedge in the corresponding triangle face.
Parameters
----------
e : edge
The edge whose opposite half-angle tangent will be calculated.
Returns
-------
tan_alpha_kij_half : double
The tangent of half the angle that lies opposite edge `e` in the
corresponding triangle.
Raises
------
ValueError
If the edge is a boundary edge, i.e. does not belong to a face.
If the triangle is degenerated or the triangle inequality doesn't hold.
Notes
-----
* The tangents of the half angles of the triangles are computed using the
half-angle theorem. For the angle :math:`\alpha` opposite `a` in a
triangle with side length `a`, `b`, `c` it is
.. math::
\tan \frac{\alpha}{2} = \sqrt{\frac{(a-b+c)(a+b-c)}{(a+b+c)(-a+b+c)}}.
"""
if e.face is None:
raise ValueError(
"The edge does not belong to a face, thus no anglewill be calculated."
)
eij = e.len
ejk = e.nex.len
eki = e.pre.len
# if eij * ejk == 0:
# raise ValueError('The triangle ({} {} {}) is degenerated.'.format(
# e.pre.head, e.head, e.nex.head))
if abs((eij**2 + ejk**2 - eki**2) / (2 * eij * ejk)) > 1:
raise ValueError(
"The triangle inequality does not hold for the"
"triangle with edges ({}{}, {}{}, {}{}).".format(
e.pre.tail, e.pre.head, e.tail, e.head, e.nex.tail, e.nex.head
)
)
try:
tan_alpha_kij_half = math.sqrt(
((eij - ejk + eki) * (eij + ejk - eki))
/ ((eij + ejk + eki) * (-eij + ejk + eki))
)
return tan_alpha_kij_half
except ZeroDivisionError:
warn(
"The triangle ({} {} {}) is degenerated.".format(
e.pre.head, e.head, e.nex.head
)
)
[docs]def opp_angle_cot(e):
r"""
Calculates the cotangent of the angle that is lying opposite the given
halfedge in the corresponding triangle face.
Parameters
----------
e : edge
The edge whose opposite angle co-tangent will be calculated.
Returns
-------
cot_alpha : double
The cotangent of the angle alpha opposite the given edge e.
Notes
-----
* The cotangent is calculated using the half-angle tangents
.. math:: \cot(\alpha) = \frac{1-\tan^2(\alpha /2)}{2\tan(\alpha /2)}.
* These cotangents can be used to form the so called cotan weight
:math:`w_{ij} = \cot(\alpha^k_{ij}) + \cot(\alpha^l_{ij})` for an edge
eij with the adjacent triangles `ijk` and `ilj`.
An edge is Delaunay, iff its cotan weight is non-negative.
"""
cot_alpha = (1 - (opp_halfangle_tan(e)) ** 2) / (2 * (opp_halfangle_tan(e)))
return cot_alpha
[docs]def euclidean_cotan_weight(e):
r"""
Calculates the Eulcidean cotan weight for an edge of a triangulation.
Parameters
----------
e : edge
The edge whose cotan weight shall be calculated.
Returns
-------
wij : double
The cotan weight of the given edge.
Notes
-----
* The Euclidean cotan weight of an edge `eij` with adjacent triangles `ijk`
and `ilj` is :math:`w_{ij} = \cot(\alpha^k_{ij}) + \cot(\alpha^l_{ij})`.
* For an edge `eij` opposite to the boundary in triangle `ijk` we set
:math:`w_{ij} = \cot(\alpha^k_{ij})`.
* The cotangents in the cotan weights are computed using the half angle
tangents of the angles opposite the edge that is to be checked.
"""
if e.opp.face is None:
return opp_angle_cot(e)
else:
wij = opp_angle_cot(e) + opp_angle_cot(e.opp)
return wij
# Version using the angle criterion instead the halfangle-tangent:
#
# alpha = calculate_angle(e.nex)
# if e.opp.face is None:
# return alpha <= math.pi
# else:
# beta = calculate_angle(e.opp.nex)
# return alpha + beta <= math.pi
[docs]def evaluate_cotan_weight(e):
return euclidean_cotan_weight(e) >= 0.0
[docs]def is_delaunay(edge_or_surface, evaluation_fct=evaluate_cotan_weight):
r"""
Checks if an edge adjacent to two triangles or a triangulated surface
is Delaunay.
Parameters
----------
edge_or_surface : edge or surface
The edge or surface that is checked for Delaunay property.
evaluation_fct : function(edge) -> bool
A function that takes an edge as argument and determines how being
delaunay is evaluated.
By default the Euclidean cotan weight is used to evaluate the delaunay
property.
Returns
-------
bool
Raises
------
ValueError
If the input is a surface that is not a triangulation.
If the input is an edge whose face is not a triangle.
Notes
-----
* By default the Delaunay property is checked by using the Euclidean
cotan weights :math:`w_{ij} = \cot(\alpha^k_{ij}) + \cot(\alpha^l_{ij})`
for an edge eij with adjacent triangles `ijk` and `ilj`.
For boundary edges, `is_delaunay(eij)` returns true, for an interior edge
it returns true if the cotan weight wij is nonnegative, false otherwise.
"""
if isinstance(edge_or_surface, Surface):
if not ddg.halfedge._get.is_triangulation(edge_or_surface):
raise ValueError("The surface is not a triangulation.")
for e in ddg.halfedge._get.interior_edges(edge_or_surface):
if not is_delaunay(e, evaluation_fct):
return False
return True
else:
e = edge_or_surface
if ddg.halfedge._get.count_edges_in_loop(e) != 3:
raise ValueError("The face belonging to e is not a triangle.")
return evaluation_fct(e)
[docs]def update_edgelength(e, length_function):
"""
Updates the edgelength attribute `len` of a given halfedge and its opposite
by use of a given function that determines the new length.
Parameters
----------
e : edge
The halfedge whose length attribute `len` will be updated according to
the given `length_function`.
The corresponding opposite edgelength `e.opp.len` will also be updated.
length_function : function(edge) -> float
A function that calculates a new edgelength for a given edge.
Returns
-------
None
"""
new_edgelength = length_function(e)
e.len = new_edgelength
e.opp.len = new_edgelength
[docs]def intrinsic_lengthfunction(e):
r"""
Function that calculates the new intrinsic edgelength for an edge that will
be flipped.
For calculation of the new edgelength after a flip the distance between
the new head and tail coordinates is measured.
Parameters
----------
e : edge
The edge whose edgelength after a flip will be calculated in
extrinsic manner.
Returns
-------
new_edgelength : float
The new intrinsic length of the corresponding flipped edge.
Notes
-----
The new edgelength :math:`\ell_{lk}` of the flipped edge is computed using
half-angle tangents and law of cosine:
.. math::
\tan(\frac{\alpha + \delta}{2}) =
\frac{\tan(\alpha /2) + \tan(delta /2)}
{1-\tan(\alpha /2)\tan(\delta /2)}
then
.. math::
\cos(\alpha + \delta) = \frac{1-\tan^2(\frac{\alpha + \delta}{2})}
{1+\tan^2(\frac{\alpha + \delta}{2})}
and for the new length
.. math:: \ell_{lk} = \sqrt{b^2 + c^2 -2bc\cos(\alpha + \delta)}.
"""
ejk = e.nex
eki = e.pre
eil = e.opp.nex
elj = e.opp.pre
# calculate and set length in intrinsic manner using halfangle tangents
tan_halfangle = (opp_halfangle_tan(eil) + opp_halfangle_tan(eki)) / (
1 - opp_halfangle_tan(eil) * opp_halfangle_tan(eki)
)
cos_angle = (1 - tan_halfangle**2) / (1 + tan_halfangle**2)
new_edgelength = math.sqrt(
elj.len**2 + ejk.len**2 - 2 * elj.len * ejk.len * cos_angle
)
return new_edgelength
# Version with calculating the angle itself instead of halfangle tangents
# alpha = calculate_angle(e) + calculate_angle(elj)
# new_edgelength = math.sqrt(elj.len**2 + ejk.len**2
# - 2*elj.len*ejk.len* math.cos(alpha))
[docs]def extrinsic_lengthfunction(e):
"""
Function that calculates the new extrinsic edgelength for an edge that will
be flipped.
For calculation of the new edgelength after a flip the distance between
the new head and tail coordinates is measured.
Parameters
----------
e : edge
The edge whose edgelength after a flip will be calculated in
extrinsic manner.
Returns
-------
new_edgelength : float
The new extrinsic length of the corresponding flipped edge.
"""
if not hasattr(e.head, "co"):
raise ValueError(
"The edge has no coordinate thus no new edgelength "
"can be calculated in extrinsic manner."
)
ejk = e.nex
eil = e.opp.nex
new_edge_head = ejk.head
new_edge_tail = eil.head
p1 = np.array(new_edge_head.co)
p2 = np.array(new_edge_tail.co)
new_edgelength = np.linalg.norm(p1 - p2)
return new_edgelength
[docs]def flip(e, length_function=intrinsic_lengthfunction):
r"""
Flips an edge within a quadrilateral of two adjacent triangles to the other
diagonal edge and by default updates its new edgelength in intrinsic
manner.
Parameters
----------
e : edge
The edge within a kite of two adjacent triangles that shall be flipped.
length_function : function(edge) -> float
A function that takes an edge (before it is being flipped) as argument
and calculates a new edgelength which will be assigned to the flipped
edge.
By default the Euclidean intrinsic length is used to update the length
of `e`.
Returns
-------
None
Raises
------
ValueError
If the edge is not lying between two adjacent triangles.
Notes
-----
* For an edge `e_ij` between the triangles `ijk` and `ilj` `flip(e_ij)`
will flip both halfedges `e_ij` and `e_ji` to the corresponding other
diagonal halfedges `e_lk` and `e_kl` in the quadrilateral `iljk`. And all
references are updated accordingly.
* This method performs an intrinsic edge flip by default, the flipped edge
is a geodesic segment in the corresponding surface thus the shape of the
the original embedded mesh does not change.
The new edgelength :math:`\ell_{lk}` of the flipped edge is computed
using half-angle tangents and law of cosine:
.. math::
\tan(\frac{\alpha + \delta}{2}) =
\frac{\tan(\alpha /2) + \tan(delta /2)}
{1-\tan(\alpha /2)\tan(\delta /2)}
then
.. math::
\cos(\alpha + \delta) = \frac{1-\tan^2(\frac{\alpha + \delta}{2})}
{1+\tan^2(\frac{\alpha + \delta}{2})}
and for the new length
.. math:: \ell_{lk} = \sqrt{b^2 + c^2 -2bc\cos(\alpha + \delta)}.
"""
if ddg.halfedge._get.is_boundary_edge(e) or ddg.halfedge._get.is_boundary_edge(
e.opp
):
raise ValueError(
"The edge lies opposite to a boundary edge and thuscannot be flipped ."
)
if ddg.halfedge._get.is_boundary_edge(e.opp):
raise ValueError("The edge is a boundary edge and thus cannot be flipped .")
if ddg.halfedge._get.count_edges_in_loop(e) != 3:
raise ValueError("The face belonging to e is not a triangle.")
if ddg.halfedge._get.count_edges_in_loop(e.opp) != 3:
raise ValueError("The face opposite to e is not a triangle.")
# if e.face is e.opp.face:
# raise ValueError('The faces at e and e.opp are the same')
# set the edges of both faces to e and e.opp respectively
# otherwise the faces might after the flip point to edges that do not
# belong to them anymore
e.face.edge = e
e.opp.face.edge = e.opp
# give indices on the vertices and name edges along them accordingly
vi = e.opp.head
vj = e.head
vk = e.nex.head
vl = e.opp.nex.head
ejk = e.nex
eki = e.pre
eil = e.opp.nex
elj = e.opp.pre
# Update edge length attribute if existing by calculating new edge length
# before actual flip and rearrangement of references as nex and pre
if hasattr(e, "len"):
update_edgelength(e, length_function)
# assure that the vertices that were incident with e before the flip will
# point to the incoming edge that stays
vi.edge = eki
vj.edge = elj
# FLIP Action: rearrange the references of the edges, and set new heads
ddg.halfedge._modify._set_next(e, eki)
ddg.halfedge._modify._set_previous(e, eil)
ddg.halfedge._modify._set_next(e.opp, elj)
ddg.halfedge._modify._set_previous(e.opp, ejk)
e.head = vk
e.opp.head = vl
ddg.halfedge._modify._set_next(eki, eil)
ddg.halfedge._modify._set_next(elj, ejk)
# renew the face reference for both edges whose face moved away
eil.face = e.face
ejk.face = e.opp.face
[docs]def delaunay_flip_algorithm(
original_surface,
return_overlay=False,
length_function=intrinsic_lengthfunction,
evaluation_fct=evaluate_cotan_weight,
):
"""
Transforms a triangulation into a Delaunay triangulation by flipping edges.
An intrinsic Delaunay edge flipping algorithm that intrinsically flips
edges of an initial triangulation such that in the end all edges are
Delaunay and so the surface is. One can pass any desired evaluation
function for determining the delaunay property as well as a length
function for calculation of new edgelength.
The result of the intrinsic_delaunay algorithm is called the corresponding
overlay structure. This surface may be non-regular, nevertheless we call
it a triangulation and do not require strong regularity for our purpose.
Parameters
----------
original_surface : surface
The triangulated surface that shall be transferred into a Delaunay
triangulation via edge flips.
return_overlay : bool
If True, an overlay structure will be created and returned, which
stores the intrinsic Delaunay Triangulation with edge attributes
`crs` giving the crossings, i.e. a list of edge tuples of
`original_surface`.
If False, the original surface will be changed to be Delaunay and
returned. False by default.
length_function : function(edge) -> float
A function that takes an edge (before it is being flipped) as argument
and calculates a new edgelength which will be assigned to the flipped
edge.
By default the Euclidean intrinsic_lengthfunction is used to update the
length of `e`.
evaluation_fct : function(edge) -> bool
A function that takes an edge as argument and determines how being
delaunay is evaluated.
By default the Euclidean cotan weight is used to evaluate the delaunay
property.
Returns
-------
surf : surface
The resulting Delaunay triangulation after edge flipping. Depending on
the argument `return_overlay` either the original surface will be
transformed and returned or an overlay structure will be created and
returned, which stores the intrinsic Delaunay Triangulation with edge
attributes `crs` giving the crossings.
Raises
------
ValueError
If the input is a surface that is not a triangulation.
Notes
-----
* By default the Delaunay property is evaluated using the corresponding
Euclidean cotan weights. If an edge is not Delaunay it will be flipped
and its length will be updated, by default using the intrinsic new
length.
* The usual definition of triangulation, which implies strong regularity,
is to narrow for our purpose. Even if the input triangulation is strongly
regular, the result of intrinsic edge flipping may be non-regular.
Therefore we do not require triangulations to be regular.
* In case of a flip, the four neighbouring edges will be stacked for
later check.
"""
if not ddg.halfedge._get.is_triangulation(original_surface):
raise ValueError("The surface is not a triangulation.")
if return_overlay:
overlay = init_overlay_structure(original_surface)
surf = overlay
else:
surf = original_surface
q = list()
for e in ddg.halfedge._get.interior_edges(surf):
if e.opp.face is not None:
q.append(e)
flipcount = 0
while len(q) != 0:
e = q.pop()
if (not is_delaunay(e, evaluation_fct)) and (e.opp.face is not None):
logger.debug("to flip is " + str(e) + "with length e.len=" + str(e.len))
flip(e, length_function)
flipcount += 1
if hasattr(e, "type"):
e.type = "flipped"
e.opp.type = "flipped"
logger.debug("new length e.len=" + str(e.len))
if return_overlay:
update_crossings(e)
logger.debug(
"crs given in verts"
+ str(
[(edge.tail.index, edge.head.index) for edge, edgeopp in e.crs]
)
)
for e_neighbour in {e.nex, e.pre, e.opp.nex, e.opp.pre}:
if e_neighbour.opp.face is not None:
q.append(e_neighbour)
logger.info("Flipped " + str(flipcount) + " edges.")
return surf
[docs]def update_crossings(e):
"""
Updates the crossings of an edge belonging to an overlay structure after
its flip, i.e. stores all original edges that are now crossed by this
current one.
Parameters
----------
e : edge
The edge whose crossings attribute `e.crs` will be updated.
Returns
-------
None
Notes
-----
The crossings of an edge, meaning which edges of the original surface are
crossed by it, are stored in the edge attribute called `crs` as a list
of edge-tuples of the form `(edge, edge.opp)`.
This list gives the halfedges in order of their crossing from `e.tail`
to `e.head`.
"""
# Collect all halfedges that cross the left of the flipped edge
# (left side = the two edges forming of the face belonging to
# e.opp without e.opp and the vertex between them)
# by adding the crossings list of the first halfedge on left side
# then append all halfedge-tuples that go out of the vertex on that side
# then add the crossings of the second (= upper half) on the left
left = e.opp.nex.crs.copy()
if left:
first_out = left[-1][1].pre
else:
original = e.head.original_vertex.surf
first_out = list(original.edges)[e.opp.nex.opp.index]
out_edge = first_out
left.append((out_edge, out_edge.opp))
out_edge = out_edge.opp.nex
while out_edge != first_out:
left.append((out_edge, out_edge.opp))
out_edge = out_edge.opp.nex
left.extend(e.opp.nex.nex.crs)
# Collect all halfedges that cross the right side (= the two edges
# belonging to the same face as e), same procedure as on other side
right = e.nex.crs.copy()
for out_edge in ddg.halfedge._get.out_edges(e.nex.head.original_vertex):
right.append((out_edge, out_edge.opp))
right.extend(e.nex.nex.crs)
# Reverse right list to later go through it in same order as left list
# when searching for edges in lists on both sides (than it is a crossing)
right.reverse()
e.crs.clear()
for (cr_opp, cr_e) in left:
if (cr_e, cr_opp) in right:
e.crs.append((cr_opp, cr_e))
e.opp.crs = list([(cr_e, cr_opp) for (cr_opp, cr_e) in reversed(e.crs)])
[docs]def init_crossings_attr(surf):
"""
Adds a list as crossings attribute 'crs' to all edges of the surface.
In this list all crossings of an edge of the current triangulation with
edges of the original triangulation are saved. This is used for bookkeeping
of and maintaining the overlay structure during flipping.
Parameters
----------
surf : surface
The surface whose edges will be added a crossings attribute that can be
addressed via 'edge.crs'.
Returns
-------
None
"""
surf.edges.add_attribute("crs")
for e in surf.edges:
e.crs = []
[docs]def init_overlay_structure(surf):
"""
Creates an inital overlay structure for the given surface.
A copy of the given surface is made and provided with the necessary
attributes in order to function as an initial overlay structure.
Parameters
----------
surf : surface
Returns
-------
overlay : surface
The copy of the given surface with reference to the original vertices
via the attribute `original_vertex` and the edge attribute `crs` for
storing the crossings of original edges.
"""
overlay = ddg.halfedge._copy.copy(
surf, ["co"], ["len"], [], original_vertex_attr="original_vertex"
)
init_crossings_attr(overlay)
return overlay
[docs]class SurfaceWithOverlay(object):
r"""
A class for constructing a visualisable intrinsic delaunay triangulation
from a halfedge surface and its intrinsic delaunay overlay structure.
The overlay structure is the result of performing the intrinsic edge flip
algorithm on the initial surface. We call this surface a triangulation,
but it may be non-regular. Also, as the word intrinsic implies, its edges
are geodesics on the initial surface, so not straight lines in
:math:`\mathbb{R}^3` anymore.
For its visualization we need to take further steps. They are implemented
by methods of this class.
"""
def __init__(self, original_surf, overlay_surf):
"""
Construct a SurfaceWithOverlay with a bidirectional vertex map between
vertices of the original surface and its overlay.
Parameters
----------
original_surf : surf
The initial halfedge surface the edge flip algorithm
(i.e. `intrinsic_delaunay(original_surf, True)`) was performed on.
overlay_surf : surf
The intrinsic delaunay triangulation of `original_surf`.
Notes
-----
As the two surfaces `original_surf` and `overlay_surf` are related by
the intrinsic delaunay algorithm, the crossings of the overlay edges
are tuples of edges of `original_surf`.
"""
self.original = original_surf
self.overlay = overlay_surf
self._vertex_map = SurfaceWithOverlay._create_bidirectional_map(
original_surf.verts, overlay_surf.verts
)
[docs] def get_intrinsic_delaunay_surf(self):
"""
Constructs a visualisable intrinsic Delaunay triangulation from a
SurfaceWithOverlay.
Returns
-------
delaunay_surf : surf
A visualisable delaunay triangulation corresponding to a
SurfaceWithOverlay.
"""
# *** INIT THE DELAUNAY SURF (DEL) as a COPY the overlay surf (OLY)
# and establish an EDGEMAP & VERTEXMAP between their edges and vertices
# In the following we will refine this copy to make it a visualisable
# intrinsic delaunay triangulation
delaunay_surf = ddg.halfedge._copy.copy(
self.overlay, ["co"], ["len"], [], original_vertex_attr="original_vertex"
)
# Add attribute 'type' to edges and verts to portray where the nodes
# come from (e.g. if it is an original edge, flipped one or comes
# from refinement
delaunay_surf.edges.add_attribute("type")
delaunay_surf.verts.add_attribute("type", "original")
edge_map = SurfaceWithOverlay._create_bidirectional_map(
self.overlay.edges, delaunay_surf.edges
)
# Create and Copy the crossings of DEL-edges from those in OLY
delaunay_surf.edges.add_attribute("crs")
for e in delaunay_surf.edges:
e.crs = edge_map[e].crs
vertex_map = SurfaceWithOverlay._create_bidirectional_map(
self.overlay.verts, delaunay_surf.verts
)
# ********************************************************************
# *** SUBDIVISION = Refine all flipped edges
# such that for each crossing a new vertex along the edge is inserted
# store the vertices that will be added by subdivision along an edge
# in a map E -> list of verts
new_verts = dict()
# also hold the crs for all new verts in crs attribute
delaunay_surf.verts.add_attribute("crs")
visited_edges = set()
for e in list(delaunay_surf.edges):
if e in visited_edges:
continue
visited_edges.add(e)
visited_edges.add(e.opp)
# If the edge is a flipped one, thus has crossings, subdivide it
if e.crs:
(
first_edge,
last_edge,
added_vertices,
) = ddg.halfedge._modify.subdivide_edge(e, len(e.crs) + 1)
# add the corresponding crossing to each new vertex
for v, crossing in zip(added_vertices, e.crs):
v.crs = crossing
# Remember where edges and vertices came from
for v in added_vertices:
v.type = "refinement"
edge = first_edge
while edge != last_edge.opp.nex:
edge.type = "flipped"
edge.opp.type = "flipped"
edge = edge.nex
# and update the edgemap such that the respective first segment
# (DEL) points to the edge of the overlay_surf that it came
# from by subdivision
assert first_edge == e
assert edge_map[first_edge] == edge_map[e]
assert first_edge.crs == e.crs
self._connect_via_map(edge_map, first_edge, edge_map[e])
self._connect_via_map(edge_map, last_edge, edge_map[e.opp])
last_edge.crs = e.opp.crs
del edge_map[e.opp]
# Store the list of verts (DEL) that were added by subdivision
# of an edge (OLY) in the map new_verts
overlay_e = edge_map[e]
new_verts[overlay_e] = added_vertices
new_verts[overlay_e.opp] = list(reversed(added_vertices))
# *** CALCULATE COORDINATES for these crossing vertices
# via a layout process of triangle strips of the original_surf
# for each edge in the overlay:
self.calculate_crs_coordinates(overlay_e, new_verts[overlay_e])
else:
e.type = "original"
e.opp.type = "original"
overlay_e = edge_map[e]
new_verts[overlay_e] = []
new_verts[overlay_e.opp] = []
# ********************************************************************
# *** INSERT EDGES to make it a triangulation
# (by subdivision we added new verts and edges thus now a face got
# more than 3 edges and verts although it has triangle shape)
# go through triangles in OLY in order to get the vertices of the real
# triangles (those with v.crs = None) easier at hand
for f_overlay in self.overlay.faces:
# match_crs_map: (V,V) -> V to find the vertices with the same crs
match_crs_map = dict()
# split_map: V -> V will store which vertices shall be connected
split_map = dict()
f = edge_map[f_overlay.edge].face # the corresponding face in DEL
for v in ddg.halfedge._get.face_vertices(f):
# go through all vertices of f (DEL) with crossing
if v.crs:
# create a tuple holding the tail and head of the edge
# that crosses via this vertex
crs_verts = tuple(
sorted((v.crs[0].tail.index, v.crs[0].head.index))
)
# find out if this vertex should be connected to one of
# original vertices
for v_triangle in ddg.halfedge._get.face_vertices(f_overlay):
if v_triangle.index in crs_verts:
split_map[v] = vertex_map[v_triangle]
# if the crossing tuple is already in the map,
# then we now have its partner with the same crs at
# hand, so store them for splitting
try:
split_map[v] = match_crs_map[crs_verts]
except KeyError:
# KeyError meaning the crs is not yet in the map,
# then add it with the vertex we have at hand as its
# value
match_crs_map[crs_verts] = v
# After collecting all pairs of vertices with the same crs we
# want to connect, go and split the face along them to do so.
# As new faces will be added, keep track which face shall be split
face_set = set()
face_set.add(f)
for v1, v2 in split_map.items(): # a list of tuple(key, value)
# go through a copy of the face_set as it changes while
# iterating
for split_face in set(face_set):
# find out which face it is that v1 and v2 share
split_face_vertices = set(
ddg.halfedge._get.face_vertices(split_face)
)
if v1 in split_face_vertices and v2 in split_face_vertices:
new_e = ddg.halfedge._modify.split_face_at(split_face, v1, v2)
face_set.add(new_e.opp.face)
new_e.type = "refinement"
new_e.opp.type = "refinement"
continue
return delaunay_surf
# TODO tests - maybe one can use pytest fixture to reuse datastructure !!
[docs] def layout_of_crossed_faces(self, e):
"""
Finds 2D coordinates for all vertices in a strip of triangles (ORI)
that are crossed by the given edge (OLY).
Parameters
----------
e : edge
The edge for which the strip of crossed triangles (ORI) will be
layed out in the plane.
Returns
-------
uv_map : V(ORI) -> [float, float]
A map assigning each vertex of the original surface a 2D
coordinate.
"""
uv_map = dict() # stores the 2D coords found in the layout process
# start with the first triangle of the strip by placing vertex e.tail
# in the origin and one edge along the x-axis
origin = np.array([0, 0])
v0 = self._vertex_map[e.tail]
uv_map[v0] = origin
crossed_edge = e.crs[0][0]
v1 = crossed_edge.tail
v2 = crossed_edge.head
c = crossed_edge.pre.len
uv_map[v1] = np.array([c, 0])
b = crossed_edge.nex.len
alpha = calculate_angle(crossed_edge.nex)
uv_map[v2] = np.array([b * math.cos(alpha), b * math.sin(alpha)])
logger.debug(self._vertex_map[e.tail])
logger.debug([(edge.tail.index, edge.head.index) for edge, edgeopp in e.crs])
logger.debug(self._vertex_map[e.head])
logger.debug({k.index: uv_map[k] for k in uv_map.keys()})
# go through the crossed edges (ORI) and find the layout of the next
# triangle of the strip by rotating the current crossed one and
# adjusting it length (as we know all angles and edgelength)
for i in range(len(e.crs)):
crossed_edge = e.crs[i][0]
alpha = -calculate_angle(crossed_edge.opp) # - for rotating clockwise
c, s = np.cos(alpha), np.sin(alpha)
rotation_by_alpha = np.array(((c, -s), (s, c)))
v1 = uv_map[crossed_edge.tail]
logger.debug("tail: " + str(crossed_edge.tail))
logger.debug("head: " + str(crossed_edge.head))
v2 = uv_map[crossed_edge.head]
uv = (
rotation_by_alpha.dot((v2 - v1)) # get rotated vector
/ crossed_edge.len
) * crossed_edge.opp.nex.len # normalization # adjust length
# assign the coordinate to the next vertex which is that one
# opposite the crossed edge
opposite_v = crossed_edge.opp.nex.head
logger.debug("opposite: " + str(opposite_v))
uv_map[opposite_v] = v1 + uv # dont forget to base vector at v1
return uv_map
[docs] def calculate_crs_coordinates(self, e, new_verts):
"""
Calculates and assigns 3D coordinates to the new vertices of a delaunay
surface that correspond to a given edge in the overlay surface.
When we evolve a delaunay surf from an overlay surf, we subdivide
those edges that crossed original ones in order to achieve straight
segments as edges. Thereby we added new vertices to the delaunay surf
which we now calculate coordinates for.
First we lay out the strip of original triangles the overlay edge runs
through in the plane (2D). Then the coordinates of the new verts (DEL)
are calculated as intersections of this overlay edge with the crossed
ones (ORI).
The intersection will be given by its barycentric coefficients
relative to the endpoints of the intersected segments. By this we do
not only get the coordinates for the intersection in the plane but
immediately also those for the corresponding point on the surface.
Parameters
----------
e : edge (OLY - in an overlay surface)
An edge of the overlay surf for which the coordinates of its
corresponding crossings vertices shall be calculated for.
new_verts : list of verts (DEL)
The vertices in the delaunay surf that were added via subdivision
of `e` which now shall be given 3D coordinates.
Returns
-------
None
"""
# Layout of those triangles (ORI) that the given edge e (OLY) crosses.
# The resulting 2D coords are given in a uv_map
uv_map = self.layout_of_crossed_faces(e)
# Get start- and end-vertex of that strip
# start is not needed any longer as we always place it in the origin
# start_v = self._vertex_map[e.tail] # = P = [0,0]
end_v = self._vertex_map[e.head] # = Q
# Go through all edges (ORI) crossed by e (OLY) to calculate its
# intersection with e and assign the coordinate
for i in range(0, len(e.crs)):
# The vertices of the crossed edge (ORI)
p1 = e.crs[i][0].head # = A
p2 = e.crs[i][0].tail # = B
# Get the intersection of AB with PQ in form its barycentric
# coefficients with respect to A,B and P,Q, where we always placed
# the start vertex P of the strip as P=[0,0]
# This needs to be passed directly because the case P == Q, as it
# occurs for the pillow (see tests), makes the Matrix singular.
coeffs = euclidean2d.intersection_in_barycentric_coords(
uv_map[p1], uv_map[p2], [0, 0], uv_map[end_v]
)
# Assign the crossings vertex (DEL) its 3D coordinate using the
# barycentric coeffs
new_verts[i].co = np.array(
[coeffs[0] * p1i + coeffs[1] * p2i for p1i, p2i in zip(p1.co, p2.co)],
)
@staticmethod
def _connect_via_map(node_map, node1, node2):
node_map[node2] = node1
node_map[node1] = node2
@staticmethod
def _create_bidirectional_map(nodes1, nodes2):
bidirectional_map = dict()
for orig_v, overlay_v in zip(nodes1, nodes2):
SurfaceWithOverlay._connect_via_map(bidirectional_map, overlay_v, orig_v)
return bidirectional_map
# How to define property methods
# = a method that looks like an attribute to the user
# but works like a function:
#
# @property
# def vertex_map(self):
# if not self._vertex_map:
# self._init_vertex_map()
# return self._vertex_map
[docs]def add_parser(sub_parsers):
"""
Parses arguments for the Delaunay edge flip algo command line application.
Takes the user-desired input arguments of the command line and performs
the Delaunay edge flip algo on the given input triangulation with it.
Notes
-----
Performs the algorithm by calling `execute_flip_algo(args)` with the
parsed arguments.
Parameters
----------
sub_parsers : ArgumentParser
The parent argument subparser of `ddg.main()` which this parser
shall belong to.
Returns
-------
None
"""
description = (
"Compute the intrinsic (or extrinsic) Delaunay Triangulation for a "
"given triangulated mesh."
)
parser = sub_parsers.add_parser(name="delaunay", description=description)
parser.add_argument(
"file",
metavar="original_mesh",
# type=argparse.FileType('r'),
type=str,
help="File containing the original mesh, format: obj or json.",
)
parser.add_argument(
"--name",
metavar="output_name",
type=str,
default="iDT_",
help="The name of the output file.",
)
parser.add_argument(
"--path",
metavar="output_path",
type=str,
default="export",
help="Output path " + "(default: %(default)s).",
)
parser.add_argument(
"--format",
metavar="out_format",
type=str,
default="obj",
help="Output file format (obj or json) " + "(default: %(default)s).",
)
parser.add_argument(
"--extr",
action="store_true",
help="Executes the extrinsic edge flip algo.",
)
parser.add_argument(
"--refine",
action="store_true",
help="Create and return the refined surface.",
)
parser.add_argument(
"--lengthsGiven",
action="store_true",
help="Edgelength are provided and shall not be calculated from Coordinates.",
)
parser.set_defaults(func=execute_flip_algo)
# args = parser.parse_args()
# if not (args.out_format == "obj" or args.out_format == "json"):
# raise ValueError("Output format option must be either obj or json")
[docs]def execute_flip_algo(args):
"""
Performs the Delaunay edge flip algorithm on a given triangulation.
Takes a triangulation given in obj or json format as input and performs
the Delaunay edge flip algorithm according to chosen optional arguments.
The resulting flipped structure is returned in a new file of format obj
(default) or json.
Notes
-----
For performing the edge flip algorithm the given input will be
converted into halfedge datastructure. And converted back into obj or
json format after performance.
Parameters
----------
args : an object holding all parsed arguments as attributes
The input arguments handed in by a user in the command line.
args.file : obj or json file
The input triangulation that shall be transformed via the edge flip
algorithm, given in obj or json format.
Returns
-------
Writes the resulting structure to a new obj/json file.
"""
# *** create default output FILENAME ***
filename = Path(args.file).stem
suffix = Path(args.file).suffix
new_name = args.output_name
if new_name == "iDT_":
if args.extr:
new_name = "eDT_"
if args.refine:
new_name = "refined_iDT_"
args.output_name = new_name + filename
# *** CONVERSION of the input file into a halfedge surface ***
if suffix == ".json":
original_surf = io.from_json(args.file)
if not args.lengthsGiven:
ddg.halfedge._set.set_euclidean_length_attr(original_surf, attr_name="len")
elif suffix == ".obj":
original_surf = io.obj_to_surface(args.file)
if not args.lengthsGiven:
ddg.halfedge._set.set_euclidean_length_attr(original_surf, attr_name="len")
else:
raise ValueError("Format of input file must be either obj or json.")
# *** Execution of the FLIP ALGO ***
if args.extr:
delaunay_surf = delaunay_flip_algorithm(
original_surf, length_function=extrinsic_lengthfunction
)
elif args.refine:
overlay = delaunay_flip_algorithm(
original_surf, True, length_function=intrinsic_lengthfunction
)
surf_with_overlay = SurfaceWithOverlay(original_surf, overlay)
delaunay_surf = surf_with_overlay.get_intrinsic_delaunay_surf()
else:
delaunay_surf = delaunay_flip_algorithm(
original_surf, length_function=intrinsic_lengthfunction
)
# *** Writing the OUTPUT FILE ***
output_path = args.output_path
if not os.path.exists(output_path):
os.mkdir(output_path)
# JUST CREATE A JSON FILE FROM OBJ INPUT:
# original_surf = io.obj_to_surface(args.file)
# set.set_euclidean_length_attr(original_surf, attr_name='len')
# io.write_json_file(original_surf, filename=output_path +
# "/" + args.output_name,
# vertex_attrs=["co"], edge_attrs=["len"])
if args.out_format == "json":
if all((hasattr(v, "co") for v in delaunay_surf.verts)):
output = io.to_json(
delaunay_surf,
filename=output_path + "/" + args.output_name + "." + args.out_format,
vertex_attrs=["co"],
edge_attrs=["len"],
)
else:
output = io.to_json(
delaunay_surf,
filename=output_path + "/" + args.output_name + "." + args.out_format,
edge_attrs=["len"],
)
else:
output = io.surface_to_obj(
delaunay_surf,
filename=output_path + "/" + args.output_name + "." + args.out_format,
)
if __name__ == "__main__":
import doctest
doctest.testmod()