Trigonometry

Throughout, we provide type signatures for all functions using jaxtyping. We use jnp (= jax.numpy) instead of numpy, and follow JAX’s functional programming paradigm (see JAX — the sharp bits).

Safe division. To avoid division by zero in degenerate configurations (e.g. zero-area triangles), all divisors are clamped via jnp.clip(x, 1e-12).

Triangle areas and circumcenters

source

get_polygon_area


def get_polygon_area(
    pts:Float[Array, 'n_vertices 2'], # Ordered polygon vertices.
)->Float[Array, '']: # Signed area.

Signed area of a 2D simple polygon (shoelace formula).

Positive for counter-clockwise vertex ordering.

source

get_triangle_area


def get_triangle_area(
    a:Float[Array, 'dim'], b:Float[Array, 'dim'], c:Float[Array, 'dim']
)->Float[Array, '']: # Triangle area.

Unsigned area of triangle with vertices a, b, c.

get_triangle_area_from_sides


def get_triangle_area_from_sides(
    u:Float[Array, 'dim'], v:Float[Array, 'dim']
)->Float[Array, '']: # Triangle area.

Unsigned area of triangle edge vectors u, v

source

get_oriented_triangle_area


def get_oriented_triangle_area(
    a:Float[Array, 'dim'], b:Float[Array, 'dim'], c:Float[Array, 'dim']
)->Float[Array, '*']: # Scalar (dim=2) or area-weighted normal vector (dim=3).

Signed area of triangle with vertices a, b, c.

source

get_circumcenter


def get_circumcenter(
    a:Float[Array, 'dim'], b:Float[Array, 'dim'], c:Float[Array, 'dim']
)->Float[Array, 'dim']: # Circumcenter coordinates.

Circumcenter of triangle with vertices a, b, c via barycentric coordinates.

source

get_voronoi_corner_area


def get_voronoi_corner_area(
    a:Float[Array, 'dim'], b:Float[Array, 'dim'], c:Float[Array, 'dim'], zero_clip:float=1e-10
)->Float[Array, '*']:

Compute Voronoi area at corner a of triangle abc. Returns zero for a degenerate triangle.

Vector operations

source

get_cot_between_vectors


def get_cot_between_vectors(
    a:Float[Array, 'dim'], b:Float[Array, 'dim']
)->Float[Array, '']: # Cotangent value.

Cotangent of the angle between vectors a and b.

source

get_angle_between_vectors


def get_angle_between_vectors(
    a:Float[Array, 'dim'], b:Float[Array, 'dim']
)->Float[Array, '']: # Angle in radians, in $[0, \pi]$.

Unsigned angle between vectors a and b.

source

get_signed_angle_between_vectors


def get_signed_angle_between_vectors(
    a:Float[Array, '2'], b:Float[Array, '2']
)->Float[Array, '']: # Signed angle in radians, in $(-\pi, \pi]$.

Signed angle from 2D vector a to b (CCW positive).

source

get_projector


def get_projector(
    normal:Float[Array, 'dim'], # Normal vector (need not be unit length).
)->Float[Array, 'dim dim']: # Projection matrix onto the plane orthogonal to normal.

Tangent-plane projector \(P = I - \hat{n} \otimes \hat{n}\).

source

project_out_vector


def project_out_vector(
    a:Float[Array, 'dim'], # Input vector.
    b:Float[Array, 'dim'], # Normal direction to project out.
)->Float[Array, 'dim']: # Component of a orthogonal to b.

Project vector a onto the plane orthogonal to b.

source

project_on_vector


def project_on_vector(
    a:Float[Array, 'dim'], # Vector to project.
    b:Float[Array, 'dim'], # Target direction.
)->Float[Array, 'dim']: # Projection of a onto b.

Project vector a onto vector b.

key1 = jax.random.key(0)
key2 = jax.random.key(1)

a, b = jax.random.normal(key1, shape=2), jax.random.normal(key2, shape=2)
assert get_triangle_area_from_sides(a, b) == jnp.linalg.norm(jnp.cross(a, b))/2

a, b = jax.random.normal(key1, shape=3), jax.random.normal(key2, shape=3)
assert get_triangle_area_from_sides(a, b) == jnp.linalg.norm(jnp.cross(a, b))/2

source

get_tetrahedron_volume


def get_tetrahedron_volume(
    a:Float[Array, '3'], b:Float[Array, '3'], c:Float[Array, '3']
)->Float[Array, '']: # Signed volume (positive if a, b, c form a right-handed frame).

Signed volume of tetrahedron with edge vectors a, b, c from a common vertex.

Intrinsic geometry

Many mesh quantities (angles, triangle areas, …) can be computed purely intrinsically from edge lengths \(\ell_{ij}\) without explicit reference to the mesh vertex coordinates in 3d.

source

get_circumcenter_from_lengths


def get_circumcenter_from_lengths(
    la:Float[Array, ''], # Length of side opposite vertex a (i.e. edge bc).
    lb:Float[Array, ''], # Length of side opposite vertex b (i.e. edge ca).
    lc:Float[Array, ''], # Length of side opposite vertex c (i.e. edge ab).
)->Float[Array, '3']: # Normalized barycentric coordinates [lambda_a, lambda_b, lambda_c].

Circumcenter in barycentric coordinates from edge lengths.

To recover Cartesian coordinates: \(u = \lambda_a \, a + \lambda_b \, b + \lambda_c \, c\).

source

get_cotangents_from_lengths


def get_cotangents_from_lengths(
    la:Float[Array, ''], # Length of side opposite vertex a (i.e. edge bc).
    lb:Float[Array, ''], # Length of side opposite vertex b (i.e. edge ca).
    lc:Float[Array, ''], # Length of side opposite vertex c (i.e. edge ab).
)->Float[Array, '3']: # Cotangents [cot_alpha, cot_beta, cot_gamma] at vertices a, b, c.

Cotangents of interior angles from side lengths.

Uses \(\cot(\alpha) = (b^2 + c^2 - a^2) / (4 \cdot \text{area})\).

source

get_angles_from_lengths


def get_angles_from_lengths(
    la:Float[Array, ''], # Length of side opposite vertex a (i.e. edge bc).
    lb:Float[Array, ''], # Length of side opposite vertex b (i.e. edge ca).
    lc:Float[Array, ''], # Length of side opposite vertex c (i.e. edge ab).
)->Float[Array, '3']: # Angles [alpha, beta, gamma] at vertices a, b, c respectively.

Interior angles from side lengths using the law of cosines.

source

get_triangle_area_from_lengths


def get_triangle_area_from_lengths(
    la:Float[Array, ''], # Length of side opposite vertex a (i.e. edge bc).
    lb:Float[Array, ''], # Length of side opposite vertex b (i.e. edge ca).
    lc:Float[Array, ''], # Length of side opposite vertex c (i.e. edge ab).
)->Float[Array, '']: # Area of the triangle.

Triangle area from side lengths using Heron’s formula.

# Test intrinsic functions against extrinsic (vertex-based) implementations
triangles = [
    jnp.array([[0., 0.], [1., 0.], [0., 1.]]),
    jnp.array([[0., 0.], [3., 0.], [1.5, 2.]]),
    jnp.array([[0., 0., 0.], [1., 0., 0.], [0., 1., 1.]]),
    jnp.array([[1., 2.], [4., 6.], [7., 1.]]),]

for tri in triangles:
    a, b, c = tri
    la, lb, lc = jnp.linalg.norm(b - c), jnp.linalg.norm(c - a), jnp.linalg.norm(a - b)

    # Area
    assert jnp.allclose(get_triangle_area(a, b, c),
                         get_triangle_area_from_lengths(la, lb, lc), atol=1e-10)
    # Angles
    angles_ext = jnp.stack([get_angle_between_vectors(b - a, c - a),
                            get_angle_between_vectors(a - b, c - b),
                            get_angle_between_vectors(a - c, b - c)])
    assert jnp.allclose(angles_ext, get_angles_from_lengths(la, lb, lc), atol=1e-10)

    # Cotangents
    cots_ext = jnp.stack([get_cot_between_vectors(b - a, c - a),
                          get_cot_between_vectors(a - b, c - b),
                          get_cot_between_vectors(a - c, b - c)])
    assert jnp.allclose(cots_ext, get_cotangents_from_lengths(la, lb, lc), atol=1e-10)

    # Circumcenter
    cc_ext = get_circumcenter(a, b, c)
    bary = get_circumcenter_from_lengths(la, lb, lc)
    cc_int = bary[0] * a + bary[1] * b + bary[2] * c
    assert jnp.allclose(cc_ext, cc_int, atol=1e-10)

Rotation matrices and normals

source

quaternion_to_rot_mat


def quaternion_to_rot_mat(
    q:Float[Array, '4'], # Quaternion $[w, x, y, z]$ (normalized internally).
)->Float[Array, '3 3']: # Rotation matrix.

Convert unit quaternion to 3D rotation matrix.

See Quaternion rotation <https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation>_.

source

get_triangle_normal


def get_triangle_normal(
    a:Float[Array, '3'], b:Float[Array, '3'], c:Float[Array, '3']
)->Float[Array, '3']: # Unit normal vector.

Unit normal of triangle abc (right-hand rule on a -> b -> c).

source

get_perp_2d


def get_perp_2d(
    x:Float[Array, '... 2'], # Input vector(s).
)->Float[Array, '... 2']: # Perpendicular vector(s).

90-degree clockwise rotation of a 2D vector: \((x, y) \to (y, -x)\).

source

get_rot_mat


def get_rot_mat(
    theta:float, # Rotation angle in radians.
)->Float[Array, '2 2']: # Rotation matrix.

2D counter-clockwise rotation matrix for angle theta.

Barycentric coordinates

source

get_barycentric_coordinates


def get_barycentric_coordinates(
    point:Float[Array, 'dim'], # Query point.
    a:Float[Array, 'dim'], b:Float[Array, 'dim'], c:Float[Array, 'dim']
)->Float[Array, '3']: # Barycentric coordinates [lambda_a, lambda_b, lambda_c] summing to 1.

Barycentric coordinates of point w.r.t. triangle abc.

Uses the area-ratio method: each coordinate is the ratio of the sub-triangle area to the full triangle area. In 3D the point is projected onto the triangle plane (least-squares sense).

vertices = jnp.array([[0., 0.], [0., 1.], [1., 0.]])
point1 = jnp.array([0.5, 0.2])

get_barycentric_coordinates(point1, *vertices)
Array([0.3, 0.2, 0.5], dtype=float64)
vertices2 = jnp.array([[0., 0., 0.], [0., 1., 0.], [1., 0., 0.]])
point2 = jnp.array([0.5, 0.2, 0.])

get_barycentric_coordinates(point2, *vertices2)
Array([0.3, 0.2, 0.5], dtype=float64)
point3 = jnp.array([0.5, 0.2, 1.])

get_barycentric_coordinates(point3, *vertices2)
Array([0.3, 0.2, 0.5], dtype=float64)

Rodrigues rotation

source

rotate_around_axis


def rotate_around_axis(
    v:Float[Array, '3'], # Vector to rotate.
    axis:Float[Array, '3'], # Unit rotation axis.
    angle:Float[Array, ''], # Rotation angle in radians.
)->Float[Array, '3']: # Rotated vector.

Rotate 3D vector v by angle (radians) around unit axis using Rodrigues’ formula.

# test: rotating x-axis by 90° around z-axis should give y-axis
x = jnp.array([1., 0., 0.])
z = jnp.array([0., 0., 1.])
result = rotate_around_axis(x, z, jnp.array(jnp.pi/2))
assert jnp.allclose(result, jnp.array([0., 1., 0.]), atol=1e-10)

# test: rotating by 0 should return the same vector
assert jnp.allclose(rotate_around_axis(x, z, jnp.array(0.)), x, atol=1e-10)

# test: rotating around the vector itself should return it
v = jnp.array([1., 2., 3.])
axis = v / jnp.linalg.norm(v)
assert jnp.allclose(rotate_around_axis(v, axis, jnp.array(1.23)), v, atol=1e-10)