r""" Module for (local) energies. All energies in this module are defined for the
minimal required amount of vertices and assume some combinatorics.
Energies in this module are not necessarily scalar valued or positive.
"""
import numpy as np
import ddg.math.complex
# Flatness Energies
# =================
[docs]def planarity_by_diagonal_distance(p1, p2, p3, p4):
r"""Flatness energy of a quad :math:`P_1P_2P_3P_4` calculate as
the distance of its diagonals. A quad is planar, if
its diagonals intersect.
Parameters
----------
p1 : numpy.ndarray of shape (n,)
Coordinates of P1
p2 : numpy.ndarray of shape (n,)
Coordinates of P2
p3 : numpy.ndarray of shape (n,)
Coordinates of P3
p4 : numpy.ndarray of shape (n,)
Coordinates of P4
Returns
-------
energy : float
(signed) circularity energy
"""
diagonal1 = p3 - p1
diagonal2 = p4 - p2
normal = np.cross(diagonal1, diagonal2)
norm = np.linalg.norm(normal)
norm = 1 if norm == 0 else norm
nnormal = normal / norm
return np.dot(nnormal, (p2 - p1))
# Circularity Energies
# ====================
[docs]def circularity_by_angles(p1, p2, p3, p4):
r"""Circularity energy of a quad :math:`P_1P_2P_3P_4` calculated by its
angles. A quad is circular, if both pairs of opposite interior angles add up
to :math:`\pi`.
Parameters
----------
p1 : numpy.ndarray of shape (n,)
Coordinates of P1
p2 : numpy.ndarray of shape (n,)
Coordinates of P2
p3 : numpy.ndarray of shape (n,)
Coordinates of P3
p4 : numpy.ndarray of shape (n,)
Coordinates of P4
Returns
-------
energy : np.ndarray of shape (2,)
(signed) circularity energy
Notes
-----
This function does not check for coplanarity.
"""
q1 = p2 - p1
q2 = p3 - p2
q3 = p4 - p3
q4 = p1 - p4
l1 = q1 / np.linalg.norm(q1)
l2 = q2 / np.linalg.norm(q2)
l3 = q3 / np.linalg.norm(q3)
l4 = q4 / np.linalg.norm(q4)
cosa = np.arccos(l1 @ -l4)
cosb = np.arccos(l2 @ -l1)
cosc = np.arccos(l3 @ -l2)
cosd = np.arccos(l4 @ -l3)
return np.array([cosa + cosc - np.pi, cosb + cosd - np.pi])
[docs]def circularity_by_ptolemy(p1, p2, p3, p4):
r"""Circularity energy of a quad :math:`P_1P_2P_3P_4` calculated by
Ptolemy's identity. A quad is circular, if
.. math ::
D_1D_2= A_1A_2 + B_1B_2,
where :math:`D_1,D_2` are the length of the quad's diagonals
and :math:`A_1,A_2` (and :math:`B_1,B_2` respectively) are
the length of opposing sides of the quad.
Parameters
----------
p1 : numpy.ndarray of shape (n,)
Coordinates of P1
p2 : numpy.ndarray of shape (n,)
Coordinates of P2
p3 : numpy.ndarray of shape (n,)
Coordinates of P3
p4 : numpy.ndarray of shape (n,)
Coordinates of P4
Returns
-------
energy : float
(signed) circularity energy
Notes
-----
This function does not check for coplanarity.
This energy is numerically bad behaved since we subtract
roots from each other.
"""
d1 = np.linalg.norm(p3 - p1)
d2 = np.linalg.norm(p4 - p2)
e1 = np.linalg.norm(p2 - p1)
e2 = np.linalg.norm(p3 - p2)
e3 = np.linalg.norm(p4 - p3)
e4 = np.linalg.norm(p1 - p4)
return d1 * d2 - e1 * e3 - e2 * e4
[docs]def circularity_by_complex_cr(p1, p2, p3, p4):
r"""Circularity energy of a quad :math:`P_1P_2P_3P_4`
in the plane calculated by complex cross ratio. A quad is circular,
if the complex cross ratio of its vertices is real.
Parameters
----------
p1 : numpy.ndarray of shape (2,)
Coordinates of :math:`P_1`
p2 : numpy.ndarray of shape (2,)
Coordinates of :math:`P_2`
p3 : numpy.ndarray of shape (2,)
Coordinates of :math:`P_3`
p4 : numpy.ndarray of shape (2,)
Coordinates of :math:`P_4`
Returns
-------
energy : float
(signed) circularity energy
"""
v1, v2, v3, v4 = map(ddg.math.complex.to_complex, (p1, p2, p3, p4))
return ddg.math.complex.cr(v1, v2, v3, v4).imag
# Isothermicity Energy
# ====================
[docs]def isothermicity_by_dual_quads(O, A, B, C, D, E, F, G, H):
r"""
Function to compute isothermicity energy of a net f
locally around a vertex :math:`f(n1, n2)` .
If a discrete surface is isothermic, real functions :math:`\alpha_k` can
be defined as edge labels that factorize face cross-ratios:
.. math::
cr(f, f_i, f_{ij}, f_j) = \frac{\alpha_i}{\alpha_j}.
With this the relationship of non-corresponding (parallel) diagonals
of primal and Christoffel dual surface can be described as
.. math::
f^*_j - f^*_i = (\alpha_i - \alpha_j) \frac{f_{ij}-f_i}{||f_{ij}-f_i||^2}.
In particular
.. math::
\langle f^*_j - f^*_i, f_{ij}-f_i \rangle = (\alpha_i - \alpha_j).
This function returns the sum of the scalar products above iterating over
all faces adjacent to the vertex :math:`f(n1, n2)` .
The sum is zero if the surface is discrete isothermic. ::
D-----C-----B
| | |
| | |
E-----O-----A
| | |
| | |
F-----G-----H
f(ni, nj)
where O = :math:`f(n1, n2)` .
Parameters
----------
O: numpy.ndarray
A: numpy.ndarray
B: numpy.ndarray
C: numpy.ndarray
D: numpy.ndarray
E: numpy.ndarray
F: numpy.ndarray
G: numpy.ndarray
H: numpy.ndarray
Returns
-------
float
Sum of scalar products :math:`\langle f^*_j - f^*_i, f_{ij}-f_i \rangle`
around vertex :math:`f(n1, n2)` .
Notes
-----
This function assumes local :math:`\mathbb{Z}^2` combinatorics of the net
around the vertex :math:`f(n1, n2)`.
"""
(
Os,
As,
Bs,
Cs,
Ds,
Es,
Fs,
Gs,
Hs,
) = ddg.math.euclidean.christoffel_dual_vertex_star(O, A, B, C, D, E, F, G, H)
return (
np.dot(B - O, Cs - As)
+ np.dot(D - O, Es - Cs)
+ np.dot(F - O, Gs - Es)
+ np.dot(H - O, As - Gs)
)