spin_transform.py

def spin_transform_deform(mesh, rho):
	# some boilerplate to convert pycomplex mesh datastructures to GA-sparse matrix operators
	I20 = mesh.topology.incidence[2, 0]   # [F, 3] face-vertex incidence
	I21 = mesh.topology.incidence[2, 1]   # [F, 3] face-edge incidence
	I10 = mesh.topology.incidence[1, 0]   # [E, 2] edge-vertex incidence
	O10 = mesh.topology._orientation[0]   # [E, 2] edge-vertex relative orientations
	O21 = mesh.topology._orientation[1]   # [F, 3] face-edge relative orientations

	T10 = as_ga_sparse(I10, as_scalar(O10))                    # edge-vertex oriented boundary operator
	A10 = as_ga_sparse(I10, as_scalar(np.ones_like(I10) / 2))  # averages vertices over edges
	A20 = as_ga_sparse(I20, as_scalar(np.ones_like(I20) / 3))  # averages vertices over faces

	M2 = as_diag(as_scalar(mesh.primal_metric[2]))          # triangle area matrix
	M2i = as_diag(as_scalar(1 / mesh.primal_metric[2]))     # inverse triangle area matrix

	vertices = as_vector(mesh.vertices)                     # cast [Vx3] float array to [V] 1vec array
	edges = T10 * vertices                                  # edge vectors from boundary operator

	# [F x V] 1-vector dirac matrix
	# this one is a little tricky;
	# relies on the fact that I20 and I21 refer to opposing vertices and edges in each entry
	D = M2i * (0.5) * as_ga_sparse(I20, edges[I21] * O21)
	# [F x V], pseudoscalar rho matrix
	R = as_diag(as_scalar(rho).dual()) * A20
	A = D - R           # [F x V], odd-grade
	Q = ~A * M2 * A     # [V x V], even-grade, or quaternionic, symmetric-positive-definite matrix
	assert Q.subspace.equals.even_grade()
	
	# now do an eigen solve
	# start with a unit quat for each vertex
	q = context.multivector.even_grade() * np.ones(mesh.topology.n_vertices)

	q = eig_solve(Q.product(q.subspace), q)
	q = q / q.norm().mean(axis=0)   # remove average scaling

	# obtain transformed edges by sandwiching with the rotor field
	transformed_edges = (A10 * q) >> edges

	# solve for new vertex positions that match transformed edges in least-square sense
	# we seek to approximate `transformed_edges = T10 * transformed_vertices`
	A = ~T10 * T10
	b = ~T10 * transformed_edges
	transformed_vertices = linear_solve(A.product(vertices.subspace), vertices, b)
	return mesh.copy(vertices=transformed_vertices.values)