Example for purely intrinsic geometry

Hi,

I've been searching for an implementation of the heat method that would be easily usable on intrinsic geometry and your library seems exactly what I was looking for.

However, I'm wondering whether the data I'm using may not be suitable. Does the implicit geometry need to be watertight?

The real data I'm looking to apply geodesic distance computations only has local intrinsic information (notably edge lengths), and no global embedding of the vertices. I have the whole global topology, but it's made of multiple groups of local 2D embeddings, and thus I was hoping your library would do the trick (especially since it seems pretty easy to compile it to Wasm and I'm planning to use it in Javascript).

However, here's a test which I was hoping would pass but currently ends up crashing:

#include "geometrycentral/surface/surface_mesh.h"
#include "geometrycentral/utilities/mesh_data.h"
#include "geometrycentral/surface/edge_length_geometry.h"
#include "geometrycentral/surface/heat_method_distance.h"
#include <Eigen/Core>
#include <stdio.h>
#include <iostream>
#include <math.h>

using namespace geometrycentral;
using namespace geometrycentral::surface;

extern "C" {

  int main(){
    // create quad with two triangles in CCW order
    //
    //   0--3
    //   | /|
    //   |/ |
    //   1--2
    //

    const size_t F = 2;
    Eigen::MatrixX3i faces(F, 3);
    faces << 0, 1, 3,
             1, 2, 3;

    Eigen::MatrixX3d edges(F, 3);
    edges << 1, M_SQRT2, 1,
             1, 1, M_SQRT2;

    std::cout << "Faces:\n" << faces << "\n";
    std::cout << "Edges:\n" << edges << "\n";

    // create surface mesh
    std::unique_ptr<SurfaceMesh> mesh;
    mesh.reset(new SurfaceMesh(faces));
    mesh->compress();
    mesh->printStatistics();

    // create geometry data
    EdgeData<double> edgeLengths(*mesh);
    for(size_t i = 0; i < faces.rows(); ++i){
      Face f = mesh->face(i);
      Halfedge he = f.halfedge(); edgeLengths[he.edge()] = edges(i, 0);
      he = he.next(); edgeLengths[he.edge()] = edges(i, 1);
      he = he.next(); edgeLengths[he.edge()] = edges(i, 2);
    }
    std::cout << "Edges:\n" << edgeLengths.raw() << "\n";
    for(Edge e : mesh->edges()){
      printf("Edge #%zu = %g\n", e.getIndex(), edgeLengths[e]);
    }
    std::unique_ptr<EdgeLengthGeometry> geometry;
    geometry.reset(new EdgeLengthGeometry(*mesh, edgeLengths));

    std::cout << "Precomputation\n";

    // create heat method distance solver (precomputation happens here)
    std::unique_ptr<HeatMethodDistanceSolver> heatSolver;
    heatSolver.reset(new HeatMethodDistanceSolver(*geometry));

    std::cout << "Solve\n";

    // solve for distance from first vertex
    const Vertex v = mesh->vertex(0);
    VertexData<double> distToSource = heatSolver->computeDistance(v);
    std::cout << "Distances:\n" << distToSource.raw() << "\n";

    return 0;
  }

}

My output (through node using Wasm) shows this:

Faces:
0 1 3
1 2 3
Edges:
      1 1.41421       1
      1       1 1.41421
Halfedge mesh with: 
    # verts =  4
    # edges =  5
    # faces =  2
    # halfedges =  6  (6 interior, 0 exterior)
      and 0 boundary components. 
Edges:
      1
1.41421
      1
      1
      1
      0
      0
      0
Edge #0 = 1
Edge #1 = 1.41421
Edge #2 = 1
Edge #3 = 1
Edge #4 = 1
Precomputation
Solve
exception thrown: 5310560

Thus the face and edge data seems appropriate (except possibly for the edges data having a few zeros for additional edges, but I guess that's due to some memory allocation strategy). The precomputation seems also fine (initially, I erroneously used M_SQRT1_2, which leads to issues during the precomputation already, but I expect that's due to the math behind since such shorter edge may lead to trouble during factorization).

My current expectation is that the distances on the square are going to be 1 from the source, except for the diagonal that ends up being again 1.41421. Am I making wrong assumptions on the data structure or expectations as to what the implicit geometry requires?

Thanks for any clarifications.

Hi! If you're looking for intrinsic heat method, this is definitely the place to be. I'm interested to hear what your use-case for purely intrinsic geometry is!

Anyway, I have two solutions for you:

First: running your code in C++, it gives the error:

std::runtime_error: ERROR: Tangent spaces not implemented for general SurfaceMesh, use ManifoldSurfaceMesh

Swapping your code to load a ManifoldSurfaceMesh rather than just a SurfaceMesh gets it to run as written. All your geometry seems reasonable, this seems to be just a datatype problem.

#include "geometrycentral/surface/manifold_surface_mesh.h"
std::unique_ptr<ManifoldSurfaceMesh> mesh;
mesh.reset(new ManifoldSurfaceMesh(faces));

Second: This error seems unreasonable, the heat method should work on nonmanifold meshes. I went ahead and made a few tweaks to allow this heat method implementation to work on nonmanifold inputs (977641442f4d43d10548cd8bdfa0d375d17e5cec); if you update to master your original code should now work as written. Also, that commit added a builtin option for more robust internal computation which might be useful to you, check it out in the docs.

So directly, the answer to your question:

Does the implicit geometry need to be watertight?

is that yesterday, our Heat Method for Distance implementation did require manifoldness, but this was due to silly code reasons. As of today, it does not require manifoldness, and furthermore has a robustification option that will improve accuracy on poorly-triangulated (possibly nonmanifold) inputs.

Let me know if you have any further questions, or this doesn't seem to work as expected!

Thanks for this very quick update! I will try the updated master, and especially look at whether it works with the geometric cases I have in practice.

The context for the fully implicit geometry is an in-browser sketch editor for describing implicit geometry that represents garments for industrial weft knitting. My sketches are combinations of edge-linked 2D polygons (because it's typically easier to draw 2D sketches than full 3D geometries). Each polygon gets its own separate coordinate system and the linking is purely topological (through edges between distinct polygons).

I got it to work with my data structure on real data. Thanks for the fixes! I guess that this closes the issue, unless there are more subtle issues. Though I have some follow-up questions due to biases I'm visually seeing in the geodesic distance distribution.

Quads vs triangles vs mixed

My original on-the-fly meshes are mostly made of a mix of quads (for the interior of polygons) and triangles (at the boundaries). Currently, I assumed I should be working with triangles only (since there are some warnings about some algorithms requiring triangles), and thus I split my quads into triangles.

Is that necessary for the heat method, or can I use mixed-triangular-quad meshes?

Impact of triangle distribution

I seem to end up with a geodesic distance that is biased (likely by the triangle directions when splitting the quads). Is that expected? Does the solution depend a lot on the shape distribution of the triangles?

I tried using the robust mode but that didn't improve the bias I am currently seeing. I will try to split my quads not into two triangles, but four triangles, which hopefully is the (main) source of the bias.

Impact of boundaries

When I run the base program with the simple quads split into two triangles, the distance I get is

d[v#0] = 0
d[v#1] = 0.7071067720614872
d[v#2] = 1.414213553248034
d[v#3] = 0.7071067720614868

The diagonal showing up as 1.41 is perfect. Couldn't be better.

However, the two corners showing up as 0.7 seems a bit odd. Is that purely due to the triangulation or are the boundaries involved?

If I change the base case with a unit quad divided into 4 triangles, then I get the following distances to the top-left (v#0):

d[v#0] = 0
d[v#1] = 1.0092490453021257
d[v#2] = 1.1955272224506772
d[v#3] = 1.0092490453021254
d[v#4] = 0.6955272224506776

The vertices are numbered as follow:

// quad with four triangles in CCW order
//
//   0---3
//   |\ /|
//   | 4 |
//   |/ \|
//   1---2
//

where v#2 is on the diagonal opposite, and v#4 is the new center vertex.

In this case, all distances are better, except for the one fully across the diagonal, that ends up lower than expected (1.19 instead of 1.41). The other ones are pretty close to expected. At least, it seems more uniform (and not skewed as I currently see in my mesh results).

Update: I tried with the updated quad splitting, and now the distribution is more uniform. The bias has seemingly disappeared, so I assume the bias was indeed an issue of triangle distribution.

Awesome! Sounds super interesting.

As you observed, the Heat Method can have some bias due to triangulation direction, since it's a PDE-based scheme. This will be exaggerated if your mesh is very coarse; the triangulation should be decently fine compared to the resolution of the geometry---for this reason you'll definitely see strange results on meshes with 4 triangles, etc.

Feel free to ask again here if you have more questions!

nmwsharp / geometry-central