For the coding portion of this assignment, you will implement various expressions for discrete curvatures and surfaces normals that you will derive in the written assignment. (However, the final expressions are given below in case you want to do the coding first.) Once implemented, you will be able to visualize these geometric quantities on a mesh. For simplicity, you may assume that the mesh has no boundary.

**Getting Started**

Please implement the following routines in `core/geometry.js`:

`angle``dihedralAngle``vertexNormalAngleWeighted``vertexNormalSphereInscribed``vertexNormalAreaWeighted``vertexNormalGaussianCurvature``vertexNormalMeanCurvature``angleDefect``totalAngleDefect``scalarMeanCurvature``circumcentricDualArea``principalCurvatures`

**Notes **

1. The dihedral angle between the normals $N_{ijk}$ and $N_{ijl}$ of two adjacent faces $ijk$ and $jil$ (respectively) is given by

$\theta_{ij} := \mathrm{atan2}( \frac{e_{ij}}{|e_{ij}|} \cdot (N_{ijk} \times N_{jil}), N_{ijk} \cdot N_{jil} )$

where $e_{ij}$ is the vector from vertex $i$ to vertex $j$.

2. The formulas for the angle weighted normal, sphere inscribed normal, area weighted normal, discrete Gaussian curvature normal and discrete mean curvature normal at vertex $i$ are

$N_i^{\phi} := \sum_{ijk \in F} \phi_i^{jk} N_{ijk}$

$N_i^S := \sum_{ijk \in F} \frac{e_{ij} \times e_{ik}}{|e_{ij}|^2 |e_{ik}|^2}$

$N_i^A := \sum_{ijk \in F} A_{ijk} N_{ijk}$

$KN_i = \frac{1}{2}\sum_{ij \in E} \frac{\theta_{ij}}{|e_{ij}|}e_{ij}$

$HN_i = \frac{1}{2}\sum_{ij \in E} (cot(\alpha_k^{ij}) + cot(\beta_l^{ij}))e_{ij}$

where $\phi_i^{jk}$ is the interior angle between edges $e_{ij}$ and $e_{ik}$, and $A_{ijk}$ is the area of face $ijk$. Note that sums are taken only over elements (edges or faces) containing vertex \(i\). Normalize the final value of all your normal vectors before returning them.

3. The circumcentric dual area at vertex $i$ is given by

$A_i := \frac{1}{8}\sum_{ijk \in F} |e_{ik}|^2 cot(\alpha_j^{ki}) + |e_{ij}|^2 cot(\beta_k^{ij})$

4. The discrete scalar Gaussian curvature (also known as angle defect) and discrete scalar mean curvature at vertex $i$ are given by

$K_i := 2\pi – \sum_{ijk \in F} \phi_i^{jk}$

$H_i := \frac{1}{2}\sum_{ij \in E} \theta_{ij} |e_{ij}|$

Note that these quantities are discrete 2-forms, i.e., they represent the total Gaussian and mean curvature integrated over a neighborhood of a vertex. They can be converted to pointwise quantities (i.e., discrete 0-forms at vertices) by dividing them by the circumcentric dual area of the vertex (i.e., by applying the discrete Hodge star).

5. You are required to derive expressions for the principal curvatures $\kappa_1$ and $\kappa_2$ in exercise 4 of the written assignment. Your implementation of `principalCurvatures` should return the (pointwise) minimum and maximum principal curvature values at a vertex (in that order).

**Submission Instructions**

Please rename your `geometry.js` file to `geometry.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 ` DDG17A2`.

Perhaps I am misinterpreting things, but I believe there are some typos in the formula for dihedral angle. At a glance, it seems like it should be sufficient to take the ratio of the norms of the cross and inner product (of the normals), probably being careful with signs/orientation. I suppose dotting with e_{i,j} helps us take care of this (assuming it is the unit vector), but it seems like the “,” should be “/”, and I don’t see where the factor of 2 comes from. Can somebody else make sense of it as written?

atan2 is the name of the arctangent function in most programming languages. Take a look at the docs here: https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/Math/atan2

Yes indeed: atan2 is a version of the arc-tangent function (i.e., inverse tangent) that takes sign into account. Wikipedia also has a thorough explanation.

Thanks for the clarification, I guess my lack of experience with various languages is showing. I’m still not getting vertexNormalGaussCurvature to work, although I’m not sure if it’s because my implementation of dihedralAngle is broken or not. Based on a comparison to an analogous formula on Wikipedia (https://en.wikipedia.org/wiki/Dihedral_angle), it seems like $e_{i,j}$ should be a unit vector (and this seems to be the most reasonable thing to do because the other factors, $N_{ijk} \times N_{ijl}, N_{ijk} \cdot N_{ijl}$ are homogenous in the edge lengths).

I have tried implementing dihedral angle via the formula from Wikipedia as well as the one provided (with $e_{i,j}$ normalized and not, with $N_{ijl} = v_{ij} \times v_{jl}$ or $N_{ijl} = v_{ji} \times v_{il}$, since I was not sure if the change in indexing from $jil$ to denote the face to $ijl$ for the subscript in the normal was intentional). However, apparently none of these have worked. Is there a way for me to test dihedralAngle independently from vertexNormalGaussCurvature?

Hi @akwon,

Apologies for the delay—for some strange reason your comment was held for moderation (most comments just get posted immediately).

My first question would be: how do you know your implementation is “wrong?” Just by looking at it? Or by running the provided tests?

Yes, you should use the unit vector \(e_{ij}/|e_{ij}|\). I believe the blog post currently reflects this fact.

Otherwise, you just need to be very careful about orientation. I.e., which direction does \(e_{ij}\) point, and which order are you taking the cross product. Since you’re walking around the neighbors \(j\) of a vertex \(i\), you need to keep these orientations consistent, i.e., \(e_{ij}\) should always point from \(i\) to \(j\), and \(N_{ijk}\), \(N_{jil}\) should be the normals of consecutive triangles in a consistent order (always counter-clockwise, for instance).

I do not get the angle weighted normal function. Apparently it is only summing the angles around the vertex. How does it return a normal vector?

Good catch. I’ve updated the formula.

Shouldn’t the value be divided by the total (such as in an average)? Apparently we are only summing the normal vectors multiplied by the angles…

You should normalize the final value of all your normal vectors before you return them.

Additionally, should we consider the cases where there is an adjacent face missing around the vertex (a hole)?

See my reply to anobani’s question below.

when there is a missing face, do we still evaluate the angle between the half edges of the missing part, and assume there is a face that can be obtained from both half edges with a normal to be added, or we skip this angle ? for example, if we have a pyramid with a missing side face, how should the angle weighted normal be related to it ?

Your solution won’t be tested on meshes with holes or boundary, but in general, the normal is only weighted by the tip angle of the faces the vertex is incident on. The adjacentFaces() iterator in the Vertex class automatically skips over “missing faces” on meshes with holes or boundary.

I am not sure I understand the vertexNormalSphereInscribed or not.

We should loop all the edges associated with vertex i for $N_i^S , ij \in E$. What is the relationship between $e_{ij}$ and $e_{ik}$? Should it be $jk \in E$ or $ijk \in F$?

Thanks!

Question: are people actually reading the course notes where these vertex normals are explained? Or just typing the expressions into the code? If it’s the latter, you may find it *very* helpful to read Section 5.4 of the notes, which also gives alternative expressions / diagrams for these normals that should make it pretty easy to see what’s going on.

Let us know if you still have questions!

I do not see a reference to the sphere-inscribed normal in section 4.5, or a section 4.5 at all… Perhaps you meant to point us to section 5.4? I do see it mentioned in 5.4.3 at http://www.cs.cmu.edu/~kmcrane/Projects/DGPDEC/paper.pdf. However, it is not obvious that the formula given for the sphere-inscribed normal is independent of an enumeration of the neighbors of a given vertex, and none is specified. It seems most reasonable that p_j, p_{j+1} should be adjacent neighbors of p_i (in the formula from Max), in which case I think ijk is a face adjacent to i?

Whoops, yes, should be 5.4, not 4.5.

Yes, j and j+1 refer to two consecutive neighbors, i.e., j and j+1 share an edge. As you say, a less ambiguous way to write this would be

\[ N_i := \sum_{ijk \in F} \frac{e_{ij} \times e_{ik}}{|e_{ij}|^2 |e_{ik}|^2} \]

(ignoring the constant).

What is $l_{ij}$?

Length of edge ij.

I’ve updated the equation to make it consistent with the others.

Is it exactly the same as $|e_{ij}|$?

Thanks.

[updated: I have just seen that you updated the notation to keep the consistency]

yes

I am confused with 5

“5. You are required to derive expressions for the principal curvatures $\kappa_1$ and $\kappa_2$ in exercise 4 of the written assignment. Your implementation of principal curvatures should return the (pointwise) minimum and maximum principal curvature values at a vertex (in that order).”

In the case of the smooth setting we have $\kappa_1$ and $\kappa_2$ related to principal curvatures of the surface. However, in the our functions for the discrete setting, our methods are based on the summation of real values. These values are associated with dihedral angle or cotan, etc.

curvature in the discrete setting.

1) In this context, what is a pointwise minimum and maximum principal curvature for a vertex?

2) Once the vertex is surrounded by n halfedges, we can have discrete methods that iterate over n values. Should we consider each of these values as a principal curvature? What does $\kappa_{i}$ mean in this context?

The expressions for the principal curvatures at a vertex are given entirely in terms of the mean and Gaussian curvature at that vertex. These expressions hold in both the discrete and smooth setting, and for both pointwise and integrated quantities. The equations for the mean and Gaussian curvature given above in 4 are integrated values, you can convert them into pointwise values (i.e., discrete 0-forms at vertices) by dividing them by the circumcentric dual area of the vertex. You can then use these pointwise values to compute your (pointwise) principal curvatures.

Does that clarify things?

1) Well, your definition implies that a vertex will have only one pointwise value for curvature – i.e. one fraction of the integrated curvature. This makes sense, but it contradicts the idea of having a min and a max for a single vertex as it is required in the function.

2) In this case, are we treating the contribution of each vertex as a “principal curvature” in the discrete setting?

3) Is this discussion related to the experiment in the Pointwise Curvature Estimates of the textbook of the course (p.68)? If not, where can I find more info about it?

as i do understand the “the dual area of the vertex” is face specific, such as the pointwise curvature value is related to each one of the dual areas.

1) Not quite, solving exercise 4 should give you expressions for both the min and the max principal curvature.

3) Related but not the same. Section 4.4.2 describes how to estimate the Gaussian curvature by fitting a quadratic polynomial to each triangle. Under refinement of the mesh, that pointwise Gaussian curvature estimate converges to its true value.

I think there is a bug in the unit test for the

`principalCurvatures`

function. The output $(\kappa_2, \kappa_1)$ of`principalCurvatures`

should always satisfy $\kappa_2 \le \kappa_1$. However, the`solution.js`

file contains expected values of $(\kappa_2, \kappa_1)$ which violate this invariant; the values of $\kappa_1$ and $\kappa_2$ are swapped in some cases. I modified the`test.js`

file to compare $\kappa_1$ (respectively, $\kappa_2$) with the maximum (respectively, minimum) of`minPrincipalCurvatures_sol[v.index]`

and`maxPrincipalCurvatures_sol[v.index]`

and my code now passes the test.You are right. A fix has been pushed to the grading script.

Shouldn’t it be |k2|<|k1| ? Aren't we just concerned about the bend magnitude?

Just checked the test script! It’s all good 🙂

Just a quick sanity check: does e.halfedge guarantee/assume an orientation?

Each halfedge points from the vertex at its base to the vertex at the base of its next halfedge.

Just to be sure: yes, for every edge there are exactly two halfedges, and they always have opposite orientation. Moreover, the orientation of every halfedge is consistent with the orientation of the triangle it sits in. So basically you should not worry that it looks any different from the standard diagrams we’ve been drawing for halfedge meshes.

A small logistical request : can the assignment and slides posts be categorized under ‘assignments’ and ‘slides’ respectively? It’s incredibly hard to find things on the website.