Example · 02

Cubify a surface

This will illustrate how to directly interact with the minimizer API for a local-global solve.

A smooth bunny changing into a cubified bunny — The original surface and the result after alternating the local and global minimizers.

This tutorial follows YASPS’s cubify_rotation_matrix example. The first example used scene.minimizeEnergy() to solve the whole simulation. Here we construct and drive the minimizer objects ourselves.

The mesh loader, renderer, and file output are omitted. Every code block below continues the same Python file.

Understand the local–global split

We optimize two fields:

Field	Shape	Updated by
`rotation`	one `3 × 3` matrix per vertex	local minimizer
`position`	one `3 × 1` vector per vertex	global minimizer

For each minimizer, we will fix the other field (by specifically not including them in Hessian computation) and find the optimal value for the field we are currently optimizing.

Provide the surface data

First import YASPS and provide the arrays produced by your mesh loader and geometry preprocessing:

import numpy as np
from yasps import attribute, minimizer, scene

rest_vertices = np.asarray([...], dtype=np.float64).reshape((-1, 3))
triangles = np.asarray([...], dtype=np.uint32).reshape((-1, 3))
edges = np.asarray([...], dtype=np.uint32).reshape((-1, 2))
vertex_normals = np.asarray([...], dtype=np.float64).reshape((-1, 3))
vertex_areas = np.asarray([...], dtype=np.float64).reshape(-1)
cotangent_weights = np.asarray([...], dtype=np.float64).reshape(-1)

edges[e] contains the two vertex indices of edge e. cotangent_weights[e] contains its matching cotangent weight. The triangle array is only needed to render or export the result.

Create the shared attributes

Now create the scene, vertex primitive, and the two live attributes:

model = scene("cubify")
surface = model.addMesh("surface")
vertices = surface.addPrimitive("vertices", numInstances=len(rest_vertices))
position = vertices.addAttribute("position", rows=3, cols=1)
rotation = vertices.addAttribute("rotation", rows=3, cols=3)
rest_position = vertices.addConstant("rest_position", rows=3, cols=1)
normal = vertices.addConstant("normal", rows=3, cols=1)
area = vertices.addConstant("area")
mass = vertices.addConstant("mass")
last_position = vertices.addConstant("last_position", rows=3, cols=1)
velocity = vertices.addConstant("velocity", rows=3, cols=1)

Initialize them with the surface data. Each rotation starts as the identity matrix:

position.updateValue(rest_vertices)
rest_position.updateValue(rest_vertices)
normal.updateValue(vertex_normals)
area.updateValue(vertex_areas)
mass.updateValue(0.001 * vertex_areas)
last_position.updateValue(rest_vertices)
velocity.updateValue(np.zeros_like(rest_vertices))
rotation.updateValue(np.tile(np.eye(3), (len(rest_vertices), 1, 1)))

We also need constants for the energy weights and time step:

dt = model.addConstant("dt")
cubic_weight = model.addConstant("cubic_weight")
arap_weight = model.addConstant("arap_weight")
orthogonality_weight = model.addConstant("orthogonality_weight")
determinant_weight = model.addConstant("determinant_weight")
regularization_weight = model.addConstant("regularization_weight")

dt.updateValue(0.01)
cubic_weight.updateValue(10.0)
arap_weight.updateValue(1.0)
orthogonality_weight.updateValue(1e4)
determinant_weight.updateValue(1e4)
regularization_weight.updateValue(1e-12)

Gather values along each edge

The ARAP energy is evaluated once per edge, so create an edge primitive and gather the two endpoint values through its connectivity:

edge_primitive = surface.addPrimitive("edges", numInstances=len(edges))
edge_to_vertex = edge_primitive.addConnectivity("edge_to_vertex", vertices, edges, 2)
edge_weight = edge_primitive.addConstant("cotangent_weight")
edge_weight.updateValue(cotangent_weights)
edge_positions = edge_primitive.addAttribute("positions", through=edge_to_vertex, source=position)
edge_rest_positions = edge_primitive.addAttribute("rest_positions", through=edge_to_vertex, source=rest_position)
edge_rotations = edge_primitive.addAttribute("rotations", through=edge_to_vertex, source=rotation)

edge_positions, edge_rest_positions, and edge_rotations now have two rows per edge, one for each endpoint.

Write the shared ARAP energy

The ARAP term compares a deformed edge with its rest edge transformed by the relaxed matrix at each endpoint:

def edge_arap_energy(weight, positions, rotations, rest_positions, penalty):
  p0 = positions.row(0).transpose()
  p1 = positions.row(1).transpose()
  q0 = rest_positions.row(0).transpose()
  q1 = rest_positions.row(1).transpose()
  rotation0 = rotations.row(0).resize(3, 3)
  rotation1 = rotations.row(1).resize(3, 3)
  deformed_edge = p1 - p0
  rest_edge = q1 - q0
  residual0 = deformed_edge - rotation0 * rest_edge
  residual1 = deformed_edge - rotation1 * rest_edge
  return 0.5 * penalty * weight * (residual0.dot(residual0) + residual1.dot(residual1))

local_arap_energy = edge_primitive.addAttribute("local_arap_energy", computed_attribute=edge_arap_energy(edge_weight, edge_positions, edge_rotations, edge_rest_positions, arap_weight))

Write cubify energy

The cubeness term encourages a rotated normal to align with a coordinate axis:

def vertex_cubic_energy(area, rotation, normal, penalty):
  rotated_normal = rotation * normal
  axis_score = rotated_normal[0].abs() + rotated_normal[1].abs() + rotated_normal[2].abs()
  return penalty * area * axis_score

cubic_energy = vertices.addAttribute("cubic_energy", computed_attribute=vertex_cubic_energy(area, rotation, normal, cubic_weight))

Energies for rotation matrix

Because rotation is a free 3 × 3 matrix rather than an angle parameterization, add two losses that keep it close to a proper rotation:

def vertex_orthogonality_energy(area, rotation, penalty):
  residual = rotation.transpose() * rotation - attribute.identity(3)
  squared_norm = residual.row(0).dot(residual.row(0)) + residual.row(1).dot(residual.row(1)) + residual.row(2).dot(residual.row(2))
  return 0.5 * penalty * area * squared_norm

def vertex_determinant_energy(area, rotation, penalty):
  residual = rotation.determinant() - 1.0
  return 0.5 * penalty * area * residual * residual

orthogonality_energy = vertices.addAttribute("orthogonality_energy", computed_attribute=vertex_orthogonality_energy(area, rotation, orthogonality_weight))
determinant_energy = vertices.addAttribute("determinant_energy", computed_attribute=vertex_determinant_energy(area, rotation, determinant_weight))

Orthogonality makes the columns perpendicular and unit length. The determinant term favors rotations with determinant 1 rather than reflections.

Build the local minimizer

Create a minimizer and register the energies that should guide rotation:

local_minimizer = minimizer()
local_minimizer.addEnergy(cubic_energy, projection_method=2)
local_minimizer.addEnergy(local_arap_energy, targets=[rotation], projection_method=2)
local_minimizer.addEnergy(orthogonality_energy, projection_method=1)
local_minimizer.addEnergy(determinant_energy, projection_method=1)
local_minimizer.addWrt([rotation])
local_minimizer.generateHessianAndGradient()

addWrt([rotation]) defines the unknowns solved by this minimizer. Since both position and rotation are attributes instead of constants, without any specification, both of the attributes will be included in the local Hessian regardless of what is specified in addWrt() function. The addWrt() only specifies what is included in the final Hessian. This is again for the correctness of the Eigenvalue decomposition. The targets=[rotation] argument however specifies that we only want to include rotation in the local Hessian. Hence, the gradient and Hessian computation will treat position as a constant.

Write the global energies

The global minimizer uses the same ARAP expression, implicit-Euler inertia, and a tiny positional regularizer:

def vertex_inertia_energy(last_position, velocity, dt, position, mass):
  target = last_position + velocity * dt
  difference = position - target
  return 0.5 * difference.transpose() * mass * difference

def vertex_regularization_energy(position, rest_position, area, penalty):
  difference = position - rest_position
  return difference.transpose() * penalty * area * difference

Create the computed attributes. The deformation terms are scaled by dt², while inertia is not:

global_arap_energy = edge_primitive.addAttribute("global_arap_energy", computed_attribute=edge_arap_energy(edge_weight, edge_positions, edge_rotations, edge_rest_positions, arap_weight) * dt * dt)
regularization_energy = vertices.addAttribute("regularization_energy", computed_attribute=vertex_regularization_energy(position, rest_position, area, regularization_weight) * dt * dt)
inertia_energy = vertices.addAttribute("inertia_energy", computed_attribute=vertex_inertia_energy(last_position, velocity, dt, position, mass))

Build the global minimizer

The second minimizer solves only for position:

global_minimizer = minimizer()
global_minimizer.addEnergy(global_arap_energy, targets=[position], projection_method=1)
global_minimizer.addEnergy(regularization_energy, projection_method=1)
global_minimizer.addEnergy(inertia_energy, projection_method=0)
global_minimizer.addWrt([position])
global_minimizer.generateHessianAndGradient()

Again, here rotation remains live inside global_arap_energy, but the local target is position.

Take one local step

Calling computeSolution() makes the minimizer assemble its Hessian and gradient, solve its linear system, and return one direction for each attribute passed to addWrt:

rotation_direction = local_minimizer.computeSolution(tolerance=1e-6, maxIterations=20_000)[0]

The local minimizer has only one target, so index 0 is the rotation direction. As in the first example, YASPS solves Hx = g, so apply the update with a minus sign. Backtrack until the local energy decreases:

rotation_before = rotation.value.copy()
local_energy_before = local_minimizer.computeTotalEnergy()
local_step = 1.0

for _ in range(8):
  rotation.updateValue(rotation_before - local_step * rotation_direction, deepCopy=True)
  if local_minimizer.computeTotalEnergy() <= local_energy_before:
    break
  local_step *= 0.5

The cubeness energy contains absolute values, so a Newton direction can land near an L1 kink. The complete example keeps the final small update if all eight trials fail.

Take one global step

The global step has the same shape, but it updates position:

position_direction = global_minimizer.computeSolution(tolerance=1e-6, maxIterations=20_000)[0]
position_before = position.value.copy()
global_energy_before = global_minimizer.computeTotalEnergy()
global_step = 1.0

for _ in range(8):
  position.updateValue(position_before - global_step * position_direction, deepCopy=True)
  if global_minimizer.computeTotalEnergy() <= global_energy_before:
    break
  global_step *= 0.5

At this point the updated position is already visible to the local minimizer, just as the updated rotation was visible to the global minimizer.

Alternate the minimizers

For each output frame, first save the previous position. Then run the local rotation minimizer to convergence before taking one global position step. Repeat until the global direction is small:

for frame in range(NUM_FRAMES):
  last_position.updateValue(position.value, deepCopy=True)

  for _ in range(MAX_POSITION_NEWTON_ITERATIONS):
    for _ in range(MAX_ROTATION_NEWTON_ITERATIONS):
      rotation_direction = local_minimizer.computeSolution(tolerance=1e-6, maxIterations=20_000)[0]
      rotation_before = rotation.value.copy()
      local_energy_before = local_minimizer.computeTotalEnergy()
      local_step = 1.0

      for _ in range(8):
        rotation.updateValue(rotation_before - local_step * rotation_direction, deepCopy=True)
        if local_minimizer.computeTotalEnergy() <= local_energy_before:
          break
        local_step *= 0.5

      if max(abs(rotation_direction)).get() < 1e-3:
        break

    position_direction = global_minimizer.computeSolution(tolerance=1e-6, maxIterations=20_000)[0]
    position_before = position.value.copy()
    global_energy_before = global_minimizer.computeTotalEnergy()
    global_step = 1.0

    for _ in range(8):
      position.updateValue(position_before - global_step * position_direction, deepCopy=True)
      if global_minimizer.computeTotalEnergy() <= global_energy_before:
        break
      global_step *= 0.5

    if max(abs(position_direction)).get() < 1e-3:
      break

  velocity.updateValue((position.value - last_position.value) / 0.01, deepCopy=True)