Harmonic mapping

Morph one surface mesh onto another via map to reference shape (disk, sphere). Part of the pipeline for dynamic data.

In this notebook, we present some additional wrapping algorithms. So far, we have wrapped two meshes by simply projection each source mesh vertex onto the closest position on the target mesh and carrying out some on-surface smoothing. This works but is not necessarily extremely robust, and can lead to very deformed UV maps when there are large, localized deformations.

Here, we’ll implement an alternative algorithm that works by (a) mapping each mesh into a common reference in a mathematically standardized way -harmonic maps- and then (b) chaining together the map from the source mesh to the common reference and from the reference to the target mesh to get a surface-surface map.

We provide algorithms for meshes of disk, cylindrical, and spherical topology (potentially with holes). For arbitrary genus (i.e. handles, like a torus), one can take a look at hyperbolic orbifolds https://github.com/noamaig/hyperbolic_orbifolds/. Not implemented here.

Loading test data

Let’s load the test meshes from the fly midgut dataset.

mesh_initial_UV = tcmesh.ObjMesh.read_obj("datasets/movie_example/initial_uv.obj")
mesh_final_UV = tcmesh.ObjMesh.read_obj("datasets/movie_example/final_uv.obj") # this is a UV map defined for tpt 20

mesh_1 = tcmesh.read_other_formats_without_uv(f"datasets/movie_example/meshes/mesh_{str(1).zfill(2)}.ply")
mesh_2 = tcmesh.read_other_formats_without_uv(f"datasets/movie_example/meshes/mesh_{str(2).zfill(2)}.ply")

mesh_60 = tcmesh.read_other_formats_without_uv(f"datasets/movie_example/meshes/mesh_{str(60).zfill(2)}.ply")
Warning: readOBJ() ignored non-comment line 4:
  o mesh_01_cylinder_seams_uv
Warning: readOBJ() ignored non-comment line 3:
  o mesh_20_uv

Disk

Following https://libigl.github.io/libigl-python-bindings/tut-chapter4/. We first map each mesh to the unit disk using harmonic coordinates. Boundary conditions are such that the boundary loop is mapped to the unit circle isometrically (i.e. relative distances are preserved).

There is still a remaining rotation degree of freedom. This is fixed by optimizing the match of the conformal factors (i.e. the area distortion of the maps to the disk) using phase correlation.

source

map_to_disk

 map_to_disk (mesh, bnd=None, set_uvs=False)

*Map mesh to unit disk by computing harmonic UV coordinates.

The longest boundary loop of the mesh is mapped to the unit circle. Follows https://libigl.github.io/libigl-python-bindings/tut-chapter4/. Note: the disk is centered at (1/2, 1/2).

The disk rotation angle is arbitrary*

Type Default Details
mesh tcmesh.ObjMesh Mesh. Must be topologically a disk (potentially with holes),
and should be triangular.
bnd NoneType None Boundary to map to the unit circle. If None, computed automatically.
set_uvs bool False whether to set the disk coordinates as UV coordinates of the mesh.
Returns np.array 2d vertex coordinates mapping the mesh to the unit disk in [0,1]^1 
mesh = tcmesh.ObjMesh.read_obj("datasets/movie_example/plane_example.obj")

disk_coordinates = map_to_disk(mesh, set_uvs=True)
Warning: readOBJ() ignored non-comment line 4:
  o Grid
plt.triplot(*disk_coordinates.T, mesh.tris)
plt.axis("equal")
(-0.049841709467476805,
 1.049992462355594,
 -0.04994346724076043,
 1.0496419841918154)

igl.flipped_triangles(disk_coordinates, mesh.tris)
array([], shape=(0, 0), dtype=int64)

Rotational alignment

The disk can still be rotated freely. Use method from the Moebius registration paper:

“We next seek the rotation that best aligns the two parameterizations. To do so, we first sample the conformal factors of each mesh onto a regular spherical grid (Figure 6, top right). We then find the rotation that maximizes the correlation between conformal factors, via a fast spectral transform.”

source

get_rot_mat2d

 get_rot_mat2d (phi)

Get 2d rotation matrix with angle phi.

mesh_a = deepcopy(mesh)
mesh_b = deepcopy(mesh)

phi = -0.8*np.pi
mesh_b.texture_vertices = ((mesh_a.texture_vertices-np.array([0.5, 0.5]))@get_rot_mat2d(phi)+np.array([0.5, 0.5]))

# test orientation flip
flip = np.diag([-1,1])
mesh_b.texture_vertices = ((mesh_b.texture_vertices-np.array([0.5, 0.5]))@flip+np.array([0.5, 0.5]))
mesh_b.faces = [fc[::-1] for fc in mesh_b.faces]
plt.triplot(*mesh_a.texture_vertices.T, mesh_a.texture_tris)

plt.triplot(*mesh_b.texture_vertices.T, mesh_b.texture_tris)
plt.axis("equal")
(-0.049841709467476805,
 1.049992462355594,
 -0.049963925863437814,
 1.0499208937180065)

To copy over attributes -like vertex 3d position- from one mesh to another, given a common parametrization, we can use barycentric interpolation.

We compute the remaining degree of freedom - rotations of the unit disk - by aligning the conformal factors of the map to the disk.

mesh_source, mesh_target = (mesh_a, mesh_b)
q = 0.01
n_grid = 1024 
allow_flip = True

disk_uv_source = np.copy(mesh_source.texture_vertices)
disk_tris_source = mesh_source.texture_tris
disk_uv_target = np.copy(mesh_target.texture_vertices)
disk_tris_target = mesh_target.texture_tris

# compute conformal factors and clip outliers
conformal_factor_source = tcdfg.compute_per_vertex_area_distortion(disk_uv_source, disk_tris_source,
                                                                   mesh_source.vertices, mesh_source.tris)
conformal_factor_source = np.clip(conformal_factor_source, np.quantile(conformal_factor_source, q),
                                  np.quantile(conformal_factor_source, 1-q))
conformal_factor_target = tcdfg.compute_per_vertex_area_distortion(disk_uv_target, disk_tris_target,
                                                                   mesh_target.vertices, mesh_target.tris)
conformal_factor_target = np.clip(conformal_factor_target, np.quantile(conformal_factor_target, q),
                                  np.quantile(conformal_factor_target, 1-q))

# interpolate into UV square
u, v = 2*[np.linspace(0, 1, n_grid),]
UV = np.stack(np.meshgrid(u, v), axis=-1).reshape((-1, 2))
interpolated_source = tcinterp.interpolate_barycentric(UV, disk_uv_source, disk_tris_source,
                                                       conformal_factor_source, distance_threshold=np.inf)
interpolated_source = interpolated_source.reshape((n_grid, n_grid))[::-1]

interpolated_target = tcinterp.interpolate_barycentric(UV, disk_uv_target, disk_tris_target,
                                                       conformal_factor_target, distance_threshold=np.inf)
interpolated_target = interpolated_target.reshape((n_grid, n_grid))[::-1]

# compute rotational alignment
interpolated_source_polar = transform.warp_polar(interpolated_source, radius=n_grid/2-1)
interpolated_target_polar = transform.warp_polar(interpolated_target, radius=n_grid/2-1)
shifts, error, _ = registration.phase_cross_correlation(interpolated_source_polar, interpolated_target_polar,
                                                        normalization=None)
if allow_flip:
    interpolated_source_polar_flipped = interpolated_source_polar[::-1]
    shifts_flipped, error_flipped, _ = registration.phase_cross_correlation(interpolated_source_polar_flipped,
                                                    interpolated_target_polar,
                                                    normalization=None)
    if error > error_flipped:
        rot_angle = shifts_flipped[0]*np.pi/180
        rot_mat = -get_rot_mat2d(rot_angle)
        new_texture_vertices = (disk_uv_source-np.array([0.5,0.5]))@rot_mat.T
        new_texture_vertices += np.array([0.5,0.5])
    
rot_angle = shifts[0]*np.pi/180
rot_mat = get_rot_mat2d(rot_angle)
new_texture_vertices = (disk_uv_source-np.array([0.5,0.5]))@rot_mat.T
new_texture_vertices += np.array([0.5,0.5])
fig, (ax1, ax2, ax3) = plt.subplots(figsize=(10,5), ncols=3)
ax1.imshow(interpolated_source_polar)
ax2.imshow(interpolated_source_polar_flipped)
ax3.imshow(interpolated_target_polar)

print()

source

rotational_align_disk

 rotational_align_disk (mesh_source, mesh_target, disk_uv_source=None,
                        disk_uv_target=None, allow_flip=True, q=0.01,
                        n_grid=1024)

*Rotationally align two UV maps to the disk by the conformal factor.

Computes aligned UV coordinates. Assumes that the UV coordinates are in [0,1]^2. This works by computing the conformal factor (how much triangle size changes as it is mapped to the plane, and finding the optimal rotation to align the conformal factors via phase correlation.*

Type Default Details
mesh_source tcmesh.ObjMesh Mesh. Must be topologically a disk (potentially with holes),
and should be triangular.
mesh_target tcmesh.ObjMesh Mesh. Must be topologically a disk (potentially with holes),
and should be triangular.
disk_uv_source NoneType None Disk coordinates for each vertex in the source mesh. Optional.
If None, the UV coordinates of mesh_source are used.
disk_uv_target NoneType None Disk coordinates for each vertex in the target mesh. Optional.
If None, the UV coordinates of disk_uv_target are used.
allow_flip bool True Whether to allow flips (improper rotations). If a flip
occurs, np.linalg.det(rot_mat) < 0
q float 0.01 Conformal factors are clipped at this quantile to avoid outliers.
n_grid int 1024 Grid for interpolation of conformal factor during alignment
Returns np.array, np.array, float new_texture_vertices_mesh_source : np.array
Rotationally aligned texture vertices
rot_mat : np.array of shape (2,2)
Rotation matrix
overlap : float
Overlap between aligned conformal factors. 1 = perfect alignment.		 
mesh_a = deepcopy(mesh)
mesh_b = deepcopy(mesh)

phi = -0.8*np.pi
mesh_b.texture_vertices = ((mesh_a.texture_vertices-np.array([0.5, 0.5]))@get_rot_mat2d(phi)+np.array([0.5, 0.5]))

flip = np.diag([-1,1])
mesh_b.texture_vertices = ((mesh_b.texture_vertices-np.array([0.5, 0.5]))@flip+np.array([0.5, 0.5]))
mesh_b.faces = [fc[::-1] for fc in mesh_b.faces]
new_vertices, rot_mat, overlap = rotational_align_disk(mesh_a, mesh_b, allow_flip=True)

overlap
0.9995530944863507
fig, (ax1, ax2, ax3) = plt.subplots(figsize=(10,5), ncols=3)
ax1.triplot(*mesh_a.texture_vertices.T, mesh_a.texture_tris)
ax2.triplot(*mesh_b.texture_vertices.T, mesh_b.texture_tris)
ax3.triplot(*new_vertices.T, mesh_b.texture_tris)
ax3.triplot(*mesh_b.texture_vertices.T, mesh_b.texture_tris)


for ax in [ax1, ax2, ax3]:
    ax.axis("equal")

Mesh registration

source

wrap_coords_via_disk

 wrap_coords_via_disk (mesh_source, mesh_target, disk_uv_source=None,
                       disk_uv_target=None, align=True, q=0.01,
                       n_grid=1024)

*Map 3d coordinates of source mesh to target mesh via a disk parametrization.

Disk parametrization can be provided or computed on the fly via harmonic coordinates. If desired, the two disks are also rotationally aligned.*

Type Default Details
mesh_source tcmesh.ObjMesh Mesh. Must be topologically a disk (potentially with holes),
and should be triangular.
mesh_target tcmesh.ObjMesh Mesh. Must be topologically a disk (potentially with holes),
and should be triangular.
disk_uv_source NoneType None Disk coordinates for each vertex in the source mesh. Optional.
If None, computed via map_to_disk.
disk_uv_target NoneType None Disk coordinates for each vertex in target mesh. Optional.
If None, computed via map_to_disk.
align bool True Whether to rotationally align the parametrizations. If False, they are used as-is.
q float 0.01 Conformal factors are clipped at this quantile to avoid outliers.
n_grid int 1024 Grid for interpolation of conformal factor during alignment.
Higher values increase alignment precision.
Returns np.array New 3d vertex coordinates for mesh_source, lying on the surface
defined by mesh_target		 
mesh_source = tcmesh.ObjMesh.read_obj("datasets/movie_example/plane_example.obj")
mesh_target = deepcopy(mesh_source)

random_rot = stats.special_ortho_group.rvs(3)
mesh_target.vertices = mesh_target.vertices@random_rot + np.array([2, 9, 1])
Warning: readOBJ() ignored non-comment line 4:
  o Grid
mesh_source = tcmesh.ObjMesh.read_obj("datasets/movie_example/plane_example.obj")
mesh_target = deepcopy(mesh_source)

random_rot = stats.special_ortho_group.rvs(3)
mesh_target.vertices = mesh_target.vertices@random_rot + np.array([2, 9, 1])
Warning: readOBJ() ignored non-comment line 4:
  o Grid
new_coords, overlap = wrap_coords_via_disk(mesh_source, mesh_target)
np.linalg.norm(new_coords-mesh_target.vertices, axis=-1).mean()
1.3403645763230948e-15

Cylinder

Map a cylinder to the disk: fill one boundary with an extra vertex, map to the disk harmonically, use Moebius to move the extra vertex to the disk center.

source

moebius_disk

 moebius_disk (pts, b)

Compute a Moebius transformation of the disk. Moves disk origin by b. pts.shape is (…, 2)

source

complex_to_xy

 complex_to_xy (arr)

Map x+iy to (x, y). Return shape==(…,2 )

source

xy_to_complex

 xy_to_complex (arr)

Map (x,y) to x+iy. arr.shape==(…,2 )

source

polygon_centroid

 polygon_centroid (pts)

See https://en.wikipedia.org/wiki/Centroid. pts.shape is (…, 2)

source

polygon_area

 polygon_area (pts)

Polygon area via shoe-lace formula. Assuming no self-intersection. pts.shape is (…, 2)

#mesh = tcmesh.ObjMesh.read_obj("movie_example/cylinder.obj") # cylinder_clean

#mesh = tcmesh.ObjMesh.read_obj("movie_example/mesh_1_cylinder_cut.obj") # cylinder_clean
mesh = tcmesh.ObjMesh.read_obj("datasets/movie_example/mesh_1_cylinder_cut_non_centered.obj") # cylinder_clean
Warning: readOBJ() ignored non-comment line 3:
  o mesh_01
# determine the boundary 
first_boundary = igl.boundary_loop(mesh.tris)
all_boundary_edges = igl.boundary_facets(mesh.tris)
second_boundary = igl.edges_to_path(np.stack([e for e in all_boundary_edges if not e[0] in first_boundary]))[0][:-1]
# add an extra vertex and triangles
filled_tris = igl.topological_hole_fill(mesh.tris, [second_boundary])
center_second_boundary = mesh.vertices[second_boundary].mean(axis=0)
filled_vertices = np.vstack([mesh.vertices, [center_second_boundary]])
filled_vertices.shape[0], mesh.vertices.shape[0], filled_tris.shape[0]-mesh.tris.shape[0], second_boundary.shape[0]
(20002, 20001, 16, 16)
# check - looks good
mesh_filled = tcmesh.ObjMesh(vertices=filled_vertices, faces=filled_tris)
mesh_filled.write_obj("datasets/movie_example/mesh_1_cylinder_cut_non_centered_filled.obj") # cylinder_clean
# map to disk via harmonic map

bnd_uv = 1 * igl.map_vertices_to_circle(filled_vertices, first_boundary)
uv = igl.harmonic(filled_vertices, filled_tris, first_boundary, bnd_uv, 1)
# map the center of the second boundary to disk center
uv = moebius_disk(uv, -uv[-1])
uv[-1] # last point -at the center of filled cylinder opening- is mapped to center of
array([0., 0.])
plt.triplot(*uv.T, filled_tris, lw=0.25)
plt.axis("equal")
(-1.096551335368751,
 1.0954288401027532,
 -1.0993918966333553,
 1.0999552780113782)

source

map_cylinder_to_disk

 map_cylinder_to_disk (mesh, outer_boundary='longest',
                       first_boundary=None, second_boundary=None,
                       set_uvs=False, return_filled=False)

*Map cylinder mesh to unit disk by computing harmonic UV coordinates.

One boundary loop of the mesh is mapped to the circle with a diameter 1/2. The second boundary is filled by adding an extra vertex at its center, which is mapped to the disk center.

Which of the two circular boundaries is mapped to the center resp. the outer circle is set by the option outer_boundary.

The disk rotation angle is arbitrary.*

Type Default Details
mesh tcmesh.ObjMesh Mesh. Must be topologically a disk (potentially with holes),
and should be triangular.
outer_boundary str longest Boundary to map to the unit circle. If “longest”/“shortest”,
the longer/shorter one is mapped to the unit circle. If int,
the boundary containing the vertex defined by the int is used.
first_boundary NoneType None First boundary loop of cylinder. If None, computed automatically.
second_boundary NoneType None Second boundary loop of cylinder. If None, computed automatically.
set_uvs bool False whether to set the disk coordinates as UV coordinates of the mesh.
return_filled bool False Whether to return vertices and faces with filled hole.
Returns np.array, np.array, np.array uv : np.array
2d vertex coordinates mapping the mesh to the unit disk in [0,1]^1.
If filled_vertices is True, the last entry is the coordinate of the added point.
filled_vertices : np.array
3d vertices with extra vertex to fill second boundary
filled_faces : np.array
Faces with added faces to fill the second boundary.		 
mesh = tcmesh.ObjMesh.read_obj("datasets/movie_example/cylinder.obj") # cylinder_clean
uv = map_cylinder_to_disk(mesh, outer_boundary="longest")
igl.flipped_triangles(uv, mesh.tris)
array([], shape=(0, 0), dtype=int64)
plt.triplot(*uv.T, mesh.tris, lw=0.5)
plt.axis("equal")
(-0.049183316665476866,
 1.049654650031695,
 -0.04980286502601551,
 1.0499106507897504)

## let's try a more challenging example - midgut mesh with two cuts to make it a cylinder

mesh = tcmesh.ObjMesh.read_obj("datasets/movie_example/mesh_1_cylinder_cut.obj") # cylinder_clean
Warning: readOBJ() ignored non-comment line 3:
  o mesh_01
uv = map_cylinder_to_disk(mesh, outer_boundary="longest",)
uv_reverse = map_cylinder_to_disk(mesh, outer_boundary="shortest", )
plt.triplot(*uv.T, mesh.tris, lw=0.5)
plt.axis("equal")
(-0.04881346336073662,
 1.0499255590960166,
 -0.049218476987450965,
 1.0493085140972238)

plt.triplot(*uv_reverse.T, mesh.tris, lw=0.5)
plt.axis("equal")
(-0.04834689694663298,
 1.047823750534748,
 -0.04965982618548542,
 1.0499708292802348)

source

wrap_coords_via_disk_cylinder

 wrap_coords_via_disk_cylinder (mesh_source, mesh_target, q=0.01,
                                n_grid=1024)

*Map 3d coordinates of source mesh to target mesh via an annulus parametrization.

Annulus parametrization can be provided or computed on the fly via harmonic coordinates. If desired, the two disks are also rotationally aligned. Meshes must be cylindrical. Two choices exist for mapping a cylinder to the plane (depending on which boundary circle is mapped to the disk boundary resp. center). Both options are tried, and the one leading to better alignment is used.

If you already have a map to the disk, use wrap_coords_via_disk instead.*

Type Default Details
mesh_source tcmesh.ObjMesh Mesh. Must be topologically a cylinder and should be triangular.
mesh_target tcmesh.ObjMesh Mesh. Must be topologically a cylinder and should be triangular.
q float 0.01 Conformal factors are clipped at this quantile to avoid outliers.
n_grid int 1024 Grid for interpolation of conformal factor during alignment.
Higher values increase alignment precision.
Returns np.array, float new_coords : np.array
New 3d vertex coordinates for mesh_source, lying on the surface
defined by mesh_target
overlap : float
Measure of geometry overlap. 1 = perfect alignment		 
mesh_source = tcmesh.ObjMesh.read_obj("datasets/movie_example/mesh_1_cylinder_cut.obj") # cylinder_clean
mesh_target = tcmesh.ObjMesh.read_obj("datasets/movie_example/mesh_1_cylinder_cut_non_centered.obj") # cylinder_clean
Warning: readOBJ() ignored non-comment line 3:
  o mesh_01
Warning: readOBJ() ignored non-comment line 3:
  o mesh_01
new_coords, overlap = wrap_coords_via_disk_cylinder(mesh_source, mesh_target,
                                                    q=0.01, n_grid=1024)
mesh_wraped = tcmesh.ObjMesh(vertices=new_coords, faces=mesh_source.faces, texture_vertices=None)
mesh_wraped.write_obj("datasets/movie_example/mesh_1_cylinder_cut_wrapepd.obj")

Sphere

For the sphere, we follow https://www.cs.cmu.edu/~kmcrane/Projects/MobiusRegistration/paper.pdf. The Riemann mapping theorem guarantees there is a conformal map of our surface to the unit sphere. This map is unique up to Moebius transformations (inversion about a point in the unit ball, and rotations). We fix the inversions by choosing the conformal map with the least amount of area distortion, and the rotation by registering the conformal factors as functions on the unit sphere using spherical harmonics

  • Cut sphere to disk by removing a single vertex
  • Map disk to plane conformally using harmonic coordinates. Turns out the least-squares conformal map is lousy at preserving angles (at least in igl)
  • Map disk to a sphere using stereographic projection, adding the removed vertex at the north pole
  • Fix Moebius inversion by choosing a conformal map with minimal distortion using Algorithm 1 from https://www.cs.cmu.edu/~kmcrane/Projects/MobiusRegistration/paper.pdf
  • Rotational registration using spherical harmonics (see notebook 03c)
  • Interpolation of vertex positions using spherical parametrization
mesh = deepcopy(mesh_final_UV)
# this is what the map of the mesh sans north pole to the plane looks like

fig = plt.figure(figsize=(4,4))

plt.triplot(*uv.T, faces_disk, lw=0.5)

array([[ 3.79091394e-08, -3.37472892e-08],
       [ 1.14168121e-07, -8.80032609e-08],
       [ 3.26394385e-08, -2.57627128e-08],
       ...,
       [ 4.22076492e-06, -5.00112042e-06],
       [ 4.04931220e-05, -5.08178908e-06],
       [ 3.67645301e-05,  5.33854828e-05]])
np.diag(average_op.toarray())
array([0., 0., 0., ..., 0., 0., 0.])

Strereographic projection

source

stereographic_sphere_to_plane

 stereographic_sphere_to_plane (pts)

*Stereographic projection from unit sphere to plane from the north pole (0,0,1).

See https://en.wikipedia.org/wiki/Stereographic_projection. Convention: the plane is at z=0, unit sphere centered at the origin. pts should be an array of shape (…, 3)*

source

stereographic_plane_to_sphere

 stereographic_plane_to_sphere (uv)

*Stereographic projection from plane to the unit sphere from the north pole (0,0,1).

See https://en.wikipedia.org/wiki/Stereographic_projection. Convention: plane is at z=0, unit sphere centered at origin. uv should be an array of shape (…, 2)*

(np.allclose(np.linalg.norm(stereographic_plane_to_sphere(uv), axis=1), 1),
 np.allclose(stereographic_sphere_to_plane(stereographic_plane_to_sphere(uv)), uv))
(True, True)

Moebius centering algorithm

Algorithm 1 from https://www.cs.cmu.edu/~kmcrane/Projects/MobiusRegistration/paper.pdf

source

center_moebius

 center_moebius (vertices_3d, vertices_sphere, tris, n_iter_centering=10,
                 alpha=0.5)

*Apply Moeboius inversions to minimize area distortion of a map from mesh to sphere.

Implementation of Algorithm 1 from: https://www.cs.cmu.edu/~kmcrane/Projects/MobiusRegistration/paper.pdf*

Type Default Details
vertices_3d np.array of shape (n_verts, 3) 3d mesh vertices
vertices_sphere np.array of shape (n_verts, 3) Initial vertex positions on unit sphere
tris np.array of shape (n_faces, 3) and type int Faces of a triangular mesh
n_iter_centering int 10 Centering algorithm iterations.
alpha float 0.5 Learning rate. Lower values make the algorithm more stable
Returns np.array, float vertices_sphere_centered : np.array of shape (n_verts, 3)
Centered sphere coordinates
com_norm : float
Distance of sphere vertex center of mass from origin. Low values
indicate convergence of the algorithm.

Putting it all together …

source

map_to_sphere

 map_to_sphere (mesh, method='harmonic', R_max=100, n_iter_centering=20,
                alpha=0.5, set_uvs=False)

*Compute a map of mesh to the unit sphere.

First, remove one vertex (the last one), and map the resulting disk-topology mesh to the plane using least squares conformal maps. Then map the plane to the sphere using stereographic projection.

The conformal map is chosen so that area distortion is as small as possible by (a) optimizing over the scale (“radius”) of the map to the disk and (b) “centering” the map using Algorithm 1 from cs.cmu.edu/~kmcrane/Projects/MobiusRegistration/paper.pdf.

This means the map is canonical up to rotations of the sphere.*

Type Default Details
mesh tcmesh.ObjMesh Mesh. Must be topologically a sphere, and should be triangular.
method str harmonic Method for computing the map from mesh without the north pole to the plane.
Recommended: harmonic.
R_max int 100 Maximum radius to consider when computing inverse stereographic
projection. If you get weird results, try a lower value.
n_iter_centering int 20 Centering algorithm iterations. If 0, no centering is performed
alpha float 0.5 Learning rate. Lower values make the centering algorithm more stable
set_uvs bool False 
mesh = deepcopy(mesh_final_UV)
coords_sphere = map_to_sphere(mesh, n_iter_centering=10, method="harmonic")
tcmesh.ObjMesh(vertices=coords_sphere, faces=mesh.tris).write_obj("datasets/movie_example/map_to_sphere.obj")
angles_3d = igl.internal_angles(mesh.vertices, mesh.tris)
angles_sphere = igl.internal_angles(coords_sphere, mesh.tris)

np.nanmean(np.abs(angles_sphere.flatten()-angles_3d.flatten())) # map to sphere is close to conformal/angle-preserving

Rotation alignment using rotation module

For a map to the sphere, we can compute the amount of area distortion at each point on the sphere. We can then use this “conformal factor” as a signal to rotationally align two maps (from different meshes) to the sphere.

# this is what the conformal factor (area distortion) looks like plotted vs phi, theta, for an example
# (the same mesh mapped to the sphere, but arbitrarily rotated)

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4))

ax1.scatter(phi_source, theta_source, c=signal_source, s=1, vmin=-3, vmax=3)
ax2.scatter(phi_target, theta_target, c=signal_target, s=1, vmin=-3, vmax=3)

ax1.axis("equal");
ax2.axis("equal");

source

rotational_align_sphere

 rotational_align_sphere (mesh_source, mesh_target, coords_sphere_source,
                          coords_sphere_target, allow_flip=False,
                          max_l=10, n_angle=100, n_subdiv_axes=1,
                          maxfev=100)

*Rotationally align two UV maps to the sphere by the conformal factor.

Computes aligned spherical coordinates. Rotational alignment works by computing the conformal factor (how much triangle size changes as it is mapped to the sphere) and optimizing over rotations to find the one which leads to the best alignment. This works via an expansion in spherical harmonics. See tcrot.rotational_alignment*

Type Default Details
mesh_source tcmesh.ObjMesh Mesh. Must be topologically a sphere (potentially with holes),
and should be triangular.
mesh_target tcmesh.ObjMesh Mesh. Must be topologically a sphere (potentially with holes),
and should be triangular.
coords_sphere_source np.array or None Sphere coordinates for each vertex in the source mesh.
If None, the UV coordinates are interpreted as angles 2piu=phi,
2piv=theta.
coords_sphere_target np.array or None Sphere coordinates for each vertex in the source mesh.
If None, the UV coordinates are interpreted as angles 2piu=phi,
2piv=theta.
allow_flip bool False Whether to allow improper rotations with determinant -1.
max_l int 10 Maximum angular momentum. If None, the maximum value available in the input
spherical harmonics is used.
n_angle int 100 Number of trial rotation angles [0,…, 2*pi]
n_subdiv_axes int 1 Controls the number of trial rotation axes. Rotation axes are vertices of
the icosphere which can be subdivided. There will be roughly
40*4**n_subdiv_axes trial axes. This parameter has the strongest influence
on the run time.
maxfev int 100 Number of function evaluations during fine optimization.
Returns np.array, np.array, float coords_sphere_source_rotated : np.array
Rotationally aligned sphere vertices
rot_mat : np.array of shape (3,3)
Rotation matrix
overlap : float
How well the conformal factors overlap. 1 = perfect overlap.		 
# first test case: same mesh, but map to the sphere rotatated

mesh_source = deepcopy(mesh_final_UV)
mesh_target = deepcopy(mesh_final_UV)

coords_sphere_source = map_to_sphere(mesh_source)
rot_mat = stats.special_ortho_group.rvs(3)
coords_sphere_target = coords_sphere_source @ rot_mat.T
transformed, R, overlap = rotational_align_sphere(mesh_source, mesh_target,
                                                  coords_sphere_source, coords_sphere_target,
                                                  allow_flip=False, max_l=10, n_angle=100, n_subdiv_axes=1,
                                                  maxfev=100)

np.linalg.norm(R-rot_mat), overlap
(0.0004874609816821069, 0.9999997693026267)
# second test case: meshes for two frames of a movie

mesh_source = deepcopy(mesh_2) # mesh_final_UV
mesh_target = deepcopy(mesh_1) # mesh_initial_UV

coords_sphere_source = map_to_sphere(mesh_source)
coords_sphere_target = map_to_sphere(mesh_target)
transformed, R, overlap = rotational_align_sphere(mesh_source, mesh_target,
                                                  coords_sphere_source, coords_sphere_target,
                                                  allow_flip=False, max_l=10, n_angle=100, n_subdiv_axes=1,
                                                  maxfev=100)
overlap
0.9908056596783977

Interpolation

Now that we have (a) mapped both meshes to the sphere, (b) minized distortion via Moebius centering, and (c) rotationally aligned the two maps via spherical harmonics, we can interpolate the coordinates of mesh B onto mesh A.

source

wrap_coords_via_sphere

 wrap_coords_via_sphere (mesh_source, mesh_target,
                         coords_sphere_source=None,
                         coords_sphere_target=None, method='harmonic',
                         n_iter_centering=10, alpha=0.5, align=True,
                         allow_flip=False, max_l=10, n_angle=100,
                         n_subdiv_axes=1, maxfev=100)

*Map 3d coordinates of source mesh to target mesh via a sphere parametrization.

Sphere parametrizations can be provided or computed on the fly using the least-area distorting conformal map to the sphere (see map_to_sphere). If desired, the two parametrizations are also aligned with respect to 3d rotations using mesh shape, using spherical harmonics. See rotational_align_sphere for details.*

Type Default Details
mesh_source tcmesh.ObjMesh Mesh. Must be topologically a disk (potentially with holes),
and should be triangular.
mesh_target tcmesh.ObjMesh Mesh. Must be topologically a disk (potentially with holes),
and should be triangular.
coords_sphere_source NoneType None Sphere coordinates for each vertex in the source mesh. Optional.
If None, computed via map_to_sphere.
coords_sphere_target NoneType None Sphere coordinates for each vertex in the target mesh. Optional.
If None, computed via map_to_sphere.
method str harmonic Method for comuting the map from sphere without north pole to plane
n_iter_centering int 10 Centering algorithm iterations for computing map to the sphere. If 0, no centering is performed
alpha float 0.5 Learning rate for computing map to the sphere. Lower values make the centering algorithm more stable.
align bool True Whether to rotationally align the parametrizations. If False, they are used as-is.
allow_flip bool False Whether to allow improper rotations with determinant -1 for rotational alignment.
max_l int 10 Maximum angular momentum. If None, the maximum value available in the input
spherical harmonics is used.
n_angle int 100 Number of trial rotation angles [0,…, 2*pi]
n_subdiv_axes int 1 Controls the number of trial rotation axes. Rotation axes are vertices of
the icosphere which can be subdivided. There will be roughly
40*4**n_subdiv_axes trial axes. This parameter has the strongest influence
on the run time.
maxfev int 100 Number of function evaluations during fine optimization for rotational alignment.
Returns np.array, float new_coords : np.array
New 3d vertex coordinates for mesh_source, lying on the surface
defined by mesh_target
overlap : np.array
Overlap of conformal factor (area distortion) on sphere of the two meshes. Only returned if align is True.
1 indicates perfect overlap.		 
#mesh_source = deepcopy(mesh_1)
#mesh_target = deepcopy(mesh_1)

#mesh_source = deepcopy(mesh_1)
#mesh_target = deepcopy(mesh_2)

mesh_source = deepcopy(mesh_final_UV)
mesh_target = deepcopy(mesh_initial_UV)
mesh_source.vertices[-1], mesh_target.vertices[-1]
(array([ 398.323151,  570.099976, -977.260376]),
 array([457.797119, 947.645935, 547.180786]))
new_coords, overlap  = wrap_coords_via_sphere(mesh_source, mesh_target,
                                              coords_sphere_source=None, coords_sphere_target=None,
                                              n_iter_centering=20, alpha=0.5, align=True,
                                              method="harmonic",
                                              allow_flip=True, max_l=10, n_angle=100, n_subdiv_axes=1, maxfev=200)
overlap # mesh_final_UV, mesh_initial_UV 0.918 with harmonic map, 0.919 harmonic-free-boundary
0.9181743310396068
mesh_wrapped = deepcopy(mesh_source)
mesh_wrapped.vertices = new_coords
mesh_source.write_obj("datasets/movie_example/source_moebius.obj")
mesh_target.write_obj("datasets/movie_example/target_moebius.obj")
mesh_wrapped.write_obj("datasets/movie_example/wrapped_moebius.obj")

Stuff that did not work out well

Non-rigid ICP

Try it out first using trimesh, then do myself to avoid dependency. Use trimesh.registration.nricp_sumner

-> No good, takes forever and my laptop runs out of memory. Also lots of parameters to tune whose meaning I don’t know.

Laplace Beltrami descriptors

Idea - compute the eigenfunctions \(\phi_i\) of the Laplace operator \(\Delta\) on the surface, and the embed each point \(p\) on the surface as \(p\mapsto \phi_i(p)/\sqrt{\lambda_i}\).

Not suitable - the spectral shapes of the midgut example mesh become extremely degenerate.

See https://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Rustamov07.pdf, https://web.archive.org/web/20100626223753id_/http://www.cs.jhu.edu/~misha/Fall07/Notes/Rustamov07.pdf

Functional mapping

[https://people.csail.mit.edu/jsolomon/assets/fmaps.pdf] - uses Laplace Beltrami descriptors. There appears to be an existing implementation in python. Let’s try out their example notebook

It does not work that well. Ok for matching from mapping timepoints 1->2, bad for 1->20, or 20->30. I get “patchy correspondences”: the map mesh A -> mesh B is not continuous and tears up mesh A into multiple patches that are stitched together in somewhat random order on mesh B. Also kinda slow.

As Dillon pointed out, shape descriptors of Laplace-Beltrami kind fundamentally use metric information. Laplace Beltrami is equivalent to a metric. Good for isometric deformations, like joint movements, bad for non-isometric, which is often the case for us.

So in principle, this is a nice method, but not quite what we need.