Assignment 5 [Coding]: Geodesic Distance (Due 5/14)

For the coding portion of your final assignment, you will implement the heat method, which is an algorithm for computing geodesic distance on curved surfaces. All of the details you need for implementation are described in Section 3 of the paper, up through and including Section 3.2. Note that you need only be concerned with the case of triangle meshes (not polygon meshes or point clouds); pay close attention to the paragraph labeled “Choice of Timestep.”

Please implement the following routines in:

  1. projects/geodesic-distances/heat-method.js:
    1. constructor
    2. computeVectorField
    3. computeDivergence
    4. compute


  • Refer to sections 3.2 of the paper for discretizations of Algorithm 1 (page 3).
  • Recall that our Laplace matrix is positive semidefinite, which might differ from the sign convention the authors use.
  • The tests for computeVectorField and computeDivergence depend on the A and F matrices you define in your constructor. So if you fail the tests but your functions look correct, check whether you have defined the flow and laplace matrices properly.
  • Your solution should implement zero neumann boundary conditions (which are the “default behavior” of the cotan Laplacian) but feel free to tryout other Dirichlet and Neumann boundary conditions on your own.

Submission Instructions

Please rename your heat-method.js file to heat-method.txt and put it in a single zip file called This file and your solution to the written exercises should be submitted together in a single email to with the subject line DDG20A5.

8 thoughts on “Assignment 5 [Coding]: Geodesic Distance (Due 5/14)”

  1. What exactly is the flow matrix? Is it M inverse times the cotan laplace matrix? Or is that just A (in heat-method.js)?

    Additionally, on page 5 of the paper near the top of the first column, it says “heat flow can be computed by solving the symmetric positive-definite system $$(M-tL_C)u=\delta_{\gamma}$$

    what is -t? Is that a typo, and should be M inverse?

    1. The top of page 94 of the attached paper go over how to build the flow matrix as you state later in the post. Our Laplace matrix is positive semidefinite, so as we have for the past two assignments, you should negate it’s sign when implementing.
      M shouldn’t be inverted. The flow term is as is.

    2. The flow matrix isn’t directly explained in the paper, but it’s basically the same as what you did for mean curvature flow in the assignment on the Laplacian.

      In the smooth setting, heat flow is governed by the partial differential equation
      \frac{d}{dt} u(\tau) = -\Delta u(\tau).
      We can approximate the time derivative by evaluating \(u\) at both \(\tau=0\) and at some time \(\tau = t\), and then taking their difference, divided by the duration \(t\):
      \frac{u(t)-u(0)}{t} = -\Delta u(t).
      Notice that on the right-hand side, the Laplacian of the function is evaluated at time \(t\). Assume that we know the initial function \(u(0)\), which in the case of the heat method is just the initial “spike” of heat at some point \(x\), i.e., \(u(0) = \delta_x\). Then to solve for the (approximate) function at time \(t\) we can rearrange the above equation to get
      u(t)+t\Delta u(t) = u(0).
      Alternatively, we can write this same equation as
      (\mathrm{id} + t\Delta) u(t) = u(0).
      From here, we still need to discretize in space. I.e., we need to turn the Laplace operator \(\Delta\) into the cotan-Laplace matrix \(L\) used on our mesh. In this case, the previous equation can be approximated by
      (I + tL) u_t = \delta,
      where \(I\) is the identity matrix. The overall matrix \(A := I + tL\) is (I guess) what you could call the “flow matrix” or “flow operator”… but this name is really not important.

      Hope that helps!

        1. Check out page 94 of the paper on heat methods. In short it seems there is no single good choice for timestep, but heuristically something related to mean edge length works well.

Leave a Reply