For the coding portion of this 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:

`projects/geodesic-distances/heat-method.js`:`constructor``computeVectorField``computeDivergence``compute`

**Notes**

- Refer to sections 3.2 of the paper for discretizations of Algorithm 1 (page 3).
- Your solution should implement zero neumann boundary conditions 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 solution.zip. This file **and** your solution to the written exercises should be submitted together in a **single email** to Geometry.Collective@gmail.com with the subject line **DDG17A5**.

Can we use late days?

Yes, you can use your remaining late days.

I have finally passed the unit test by tweaking the pos/neg signs here and there. I finally get it working by negating the integrated divergence. May I have some clarification on this? Is it because the paper assumes that the Xj’s come from the non-negated gradients? Thanks!

The laplace matrix in the paper is negative definite. The negative sign in the integrated divergence calculation compensates for our laplace matrix being positive definite. The reference solution has a bug – it should check for this negative sign in the 3rd test, not in the 2nd for computeDivergence. Sorry about the trouble, thanks for bringing this to my notice!

No problem! Thanks for the clarification!

Speaking of which, it seems that the test script does not call the “compute” method at all, which is absurd…

The test script checks your distance values by solving the poisson equation with your integrated divergence values on the right hand side. I’m not going to assign much weight to the compute method, since it mostly calls other methods.

In the constructor, what is the “mean curvature flow” matrix? Should it be the heat flow matrix?

Yes, the operator you built for the 3rd assignment.

In the step of normalizing the vector field, should I be normalizing each vector by the norm of the sum of all the vectors over each face on the mesh? Or is it local normalization, per face?

Just per face. Think about the eikonal equation: it says that \(|\nabla \phi|=1\), i.e., the norm of the gradient of the distance function equals one

at each point. So, you want your vector field to have unit norm at each point as well.Thank you for your response.

Then, when eiknoal equation has a value in between [0,1) on the RHS, would this cause our vector field to have bigger value than 1 at any point?

Hmm… not quite sure what you mean. The eikonal equation always 1 on the right-hand side. If you try to integrate a vector field whose norm is not 1 at each point (i.e., if you try to solve \(\Delta \phi = \nabla \cdot X\) when \(|X| \ne 1\)), it’s not clear what relationship (if any) the resulting function will have to distance.

I think our constraints are starting to make more sense now, thank you for your reply!

I have just started coding. I have two problems:

1) For some reason, the constructor cannot access the class MeanCurvatureFlow, so I had to copy the code for buildFlowOperator. I am not sure if this is a good practice.

2) In computeVectorField, I traverse the faces and try to access the correspondent halfedges. However, in the first faces I get “undefined” for halfedge. I assume that this is an outer face, so I added the vector (0, 0, 0) to the dictionary. However, the solution has values such as (0.7087172948207671, 0.45343192374061964, 0.5404806070108483).

1) You’ll have to modify the index.html file to access MeanCurvatureFlow. However, you shouldn’t have to – in the Geodesic’s constructor you have access to the laplace and mass matrices. You can construct the flow operator very easily from them (just as you did in assignment 3).

2) Try accessing the halfedges in each face by using the adjacentHalfedges() iterator.

1) Yes.. that is what I did…

2) The strange thing now is that I am traversing the faces in this.geometry.mesh.faces, but when I check the type with typeof(face) it says it is a string “0”, which might explains why the .halfedge was undefined. [SOLVED] I was using “in” instead of the “of” for the loop.

I am pretty sure that my function is similar to the paper, but I get a small difference for the vector:

Solution: Vector {x: 0.7087172948207671, y: 0.45343192374061964, z: -0.5404806070108483}

Mine: Vector {x: 0.725684220863355, y: 0.47096515561476415, z: -0.5015717633466883}

The error seems to be too large to be related to the browser, but it is too small to be coincidence.

It is probably related to the time step…

Your time step should equal the squared mean edge length of the mesh. There is a function in geometry.js to compute the mean edge length.

Yes… I am working on this. Thanks for the reply.