Assignment 2 (Coding): Investigating Curvature

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:

  1. angle
  2. dihedralAngle
  3. vertexNormalAngleWeighted
  4. vertexNormalSphereInscribed
  5. vertexNormalAreaWeighted
  6. vertexNormalGaussianCurvature
  7. vertexNormalMeanCurvature
  8. angleDefect
  9. totalAngleDefect
  10. scalarMeanCurvature
  11. circumcentricDualArea
  12. 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.

35 thoughts on “Assignment 2 (Coding): Investigating Curvature”

  1. 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?

      1. 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?

        1. 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).

      1. 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 ?

        1. 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.

  2. 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!

    1. 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!

      1. 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?

        1. 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).

  3. 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?

    1. 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. 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?

        1. 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.

  4. 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.

    1. 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.

  5. 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.

Leave a Reply