Add stripe patterns algorithm

[DavidJourdan]

I've implemented the stripes algorithm in my own fork of geometry central and I thought you might be interested in adding it to the library!

The function outputs values stored in a CornerData container, I did not find a way to visualize the stripes easily using Polyscope, but we can check it works by extracting the isolines and comparing them with the original implementation.

We've been using this function for a while now and everything seems to work fine, please tell me if there's anything that needs to be changed.

[nmwsharp]

Awesome, thank you for contributing!!

I scrolled through the code and it all looks good to me. I will find some time to add docs and maybe a tutorial w/ Polyscope viz to this commit, then merge it in.

[nmwsharp]

By the way: you mentioned extracting isolines. If you have code for that which you would like to share, it would be very welcome in geometry-central! Extracting isolines is a very common need for this algorithm, but I don't believe there is code available anywhere.

[DavidJourdan]

You're right, there's some functions for extracting isolines e.g. in libigl, but since we want to extract the 0 (mod 2pi) isolines we can't directly reuse these methods. I just pushed our isolines code, Etienne @evouga also contributed to this so I have to thank him for that :)

Since you mentioned writing docs, here's a short snippet to show how to use these functions, I hope that helps:

#include "geometrycentral/surface/direction_fields.h"
#include "geometrycentral/surface/meshio.h"
#include "geometrycentral/surface/stripe_patterns.h"
using namespace geometrycentral;
using namespace surface;

// Load a mesh
std::unique_ptr<HalfedgeMesh> mesh;
std::unique_ptr<VertexPositionGeometry> geometry;
std::tie(mesh, geometry) = loadMesh(filename);

// Compute a smooth direction field on vertices
VertexData<Vector2> directions = computeSmoothestVertexDirectionField(*geometry, 2);
// Constant frequency values
VertexData<double> frequencies(*mesh, 1.0);

// compute stripe pattern values
CornerData<double> stripeValues;
FaceData<int> stripeIndices;
FaceData<int> fieldIndices;
std::tie(stripeValues, stripeIndices, fieldIndices) = computeStripePattern(*geometry, directions, frequencies);

// extract isolines from these values
std::vector<Vector3> positions;
std::vector<std::array<int, 2>> edges;
std::tie(positions, edges) = extractPolylinesFromStripePattern(*geometry, stripeValues, stripeIndices, fieldIndices);

// visualize with Polyscope
polyscope::registerSurfaceMesh("mesh", geometry->inputVertexPositions, mesh->getFaceVertexList());
polyscope::registerCurveNetwork("isolines", positions, edges);

[evouga]

Yes, there is libigl code for extracting isolines, though it has some issues that as far as I know haven't been fixed: https://github.com/libigl/libigl/issues/1245

[nmwsharp]

Hi all, thanks again for submitting this.

I went to test/merge/tutorial this, but had some trouble getting it to work on various random input meshes. I kept running in to

std::runtime_error: Isolines should only branch out on singularities

or else no isolines would be generated at all. Perhaps I'm doing something wrong. Could you share a mesh + frequency settings for which it gives good isolines?

Separately, I noticed this:

// WARNING this only works if there is no more than one isoline crossing at any edge of the geometry

Was there any big blocker for this? Is it just to avoid the difficult of matching up which isoline is which at shared edges?

[evouga]

The isoline extraction problem has three different settings:

Real Case

Given a real-valued discrete function on a mesh (i.e. one number per vertex), and a real number c, extract the isolines with isovalue c', where c' is infinitesimally larger than c (an arbitrary real number larger than c but smaller than the next representable floating-point number after c).

The isolines in this case always are 1-manifolds, cross each edge of the mesh at most once, and have no singularities, branches, or other shenanigans (the c' business guarantees this, even if one or more vertices of the mesh has value exactly equal to c). The code should be 100% robust in this setting; if not let me know as I have some older working version of the code and can investigate.

Periodic Case

Given an S^1-valued discrete function \theta on a mesh and a real number c, extract the isolines with isovalue c', where c' is infinitesimally larger than c.

Isolines are still guaranteed to cross each edge at most once, but are no longer guaranteed to be topological loops, since there can be edge dislocations at triangles where d^2\theta != 0. The code should still be robust in this setting, although, the isolines will stop at the boundaries of triangles containing singularities.

Knoppel Case

Given a real-valued discrete function (i.e. one number per vertex), a discrete one-form d\theta on the edges of the mesh encoding period jumps in \theta, and a real number c, extract the isolines with isovalue c', where c' is infinitesimally larger than c.

This is the version of the problem needed for the stripes paper; now multiple isolines can cross a given edge and the problem of connecting up isolines within triangles is not trivial (especially if the total number of crossings around the edges of a triangle is odd...)

I don't believe David has tried to implement this case fully? I know my code definitely does not handle it.

[DavidJourdan]

Hi, I think the isolines extraction function should work correctly now, but you might still be running into that std::runtime_error if your mesh is not high resolution enough.

I agree with what Etienne said, the extractIsolinesFromStripePattern function solves the periodic case but not the Knoppel case. The original stripe patterns paper describes a way to render high frequency stripes on low res meshes, but I don't think there is a simple way to translate that into an isoline extraction procedure.

However, what we can do in practice is try to detect if two or more isolines might be crossing the same triangle. To do that we can compute the difference between the values at each pair of corners, and if they differ more than 2pi that means we should probably subdivide the mesh.

double maxDifference = 0;
for(Halfedge h: mesh->halfedges())
{
  if(std::abs(stripesValues[h.corner()] - stripesValues[h.next().corner()]) > maxDistance)
    maxDifference = std::abs(stripesValues[h.corner()] - stripesValues[h.next().corner()]);
}
if(maxDifference > 2 * PI)
  // subdivide the mesh

[keenancrane]

Hi all,

If you want to handle the "Knöppel case" where many isocontours pass through the same triangle, you can basically use the algorithm described in Section 3.4.2 of Sharp et al, "Navigating Intrinsic Triangulations"—at least for triangles that do not contain singularities.

The important (and maybe surprising) observation is that if you know how many isolines cross each edge, then this already uniquely determines how the crossings must get connected up. In fact, all you need is that the number of crossings n_ij on the three edges satsify the usual triangle inequalities n_ij + n_jk ≥ n_ki, and that the sum n_ij + n_jk + n_ki is even (so that every segment coming "in" can be matched up with one going "out").

Here's a description of the algorithm where I've redacted the words "intrinsic" and "extrinsic" (from the original paper) to make it clearer / more generic:

There should be code in geometry-central (or at least the intrinsic triangulations demo) that implements this, but probably better to incorporate it directly in a generic isoline extraction routine.

[keenancrane]

The singular case is ambiguous if all you know is the number of crossings. Consider for instance this example from the Stripe Patterns paper—there are multiple valid ways to connect up the segments:

If you know which edge crossing connects to the singularity in the middle, then the algorithm becomes straightforward again: pick any corner i, and connect as many segments along the incident edges ij and jk (as in Case 2 above) as long as you are not starving the other two corners.

For a triangle containing a branch point, there will be exactly one singularity and this algorithm works. In general, however, triangles can contain multiple zeros. In this case, unfortunately, things are not as simple. You could try to devise some general isoline extraction algorithm, but I think a more straightforward approach would be to subdivide triangles until they contain at most one zero (which you can detect from the data on edges). Then apply the one-singularity extraction algorithm.

[DavidJourdan]

Thanks for the detailed explanations! I don't have a lot of time right now so it might take a while before I add these features to the isoline extraction procedure though. In the meantime the isoline extraction procedure should work fine if the mesh if high resolution enough, and it will throw an error message if there are several stripes passing through the same triangle.

[rahul-mitra13]

Awesome, thank you for contributing!!

I scrolled through the code and it all looks good to me. I will find some time to add docs and maybe a tutorial w/ Polyscope viz to this commit, then merge it in.

Hey all! It seems as though the polyscope viz was never integrated. Is there a place I can look to write up the shader code for the stripe visualization?

[DavidJourdan]

Hi, After discussing with Rahul at SIGGRAPH I decided to revive this PR and finally fix the limitation that we can't have more than one isoline passing through a given triangle. Rahul's idea was that within each non-singular triangle we can identify isolines by their associated isovalue (i.e. $2k\pi, k \in \mathbb{Z}$), which makes matching the crossings easier.

Now for singularities it becomes harder because since there is a $2k\pi$ jump, we can't directly match crossings by their associated isovalue. And as Keenan pointed out, this matching is ambiguous if all you know is the number and positions of the crossings.

What I ended up doing is matching points which are near each other and whose direction is close to perpendicular to the input directional field (since the generated stripes are perpendicular to the input field). In practice it means finding the pair of points $v, w$ which minimize $(v - w) \cdot X_i$ and iterate until all the matchings are made.

I left this as an option which is disabled by default because in practice we can often find cases where isolines are crossing each other on singularities, which shouldn't happen. The best way to avoid this might be to implement the routines from the intrinsic triangulations paper that Keenan described (I haven't found them in the intrinsic triangulations demo though).

[nmwsharp]

Hi David, very happy to hear it, it would be great to get this in! I don't have a ton of time these days, but happy to help how I can.

I like the idea of using a guiding direction field to chose from ambiguous cases. Perhaps you could put this to work in a slightly different way, but noting that the algorithm Keenan mentions above is basically an 'corner clipping' procedure---you repeatedly pick a corner and connect the two crossings along edges closest to that corner. The observation is that long as you don't choose the corner adjacent to the two edges with the fewest crossings, this will always work. The trick with your 'singular' cases is that there are will be extra crossings left over when cannot clip any more corners; these crossings then need to be connected to the central singularity. Perhaps if you use the guiding field to choose which corner to clip, you will often be left with a geometrically-desirable solution.

By the way, the related routines that Keenan mentions from the intrinsic triangulations stuff appear here, https://github.com/nmwsharp/geometry-central/blob/master/src/surface/common_subdivision.cpp#L553-L651 although it may not be easy to simply re-use them.

Also, as you're working on this feel free to let me know when it seems 'ready' on your end, or if you have any direct questions.

[DavidJourdan]

Hi, the isoline extraction and matching on singularities is mostly done on my end. Ideally it would be nice to treat the case of branch points in the input direction field separately from the other singularities, but I don't think there's any major issues with the way it is implemented right now so I think it is 'ready' :)

What I ended up doing in the end is I implemented the routines described by Keenan (ll.280-341), and to make them work in the case of singular triangles, we remove $n$ candidate crossings using the guiding field (where $n$ is the multiplicity of the singularity in absolute value) in a way that preserves the triangular inequality $n_{ij} \leq n_{jk} + n_{ki}$.

I have tested this implementation on a fork of the gc-polyscope-project-template with various meshes and it seems to work perfectly even with very high frequencies now.

[DavidJourdan]

The only last feature I would like to implement (but perhaps it's better to submit this in a separate pull request) is the face-based alternative, that is computing stripe patterns with input data residing on faces instead of vertices.

I have implemented most of the functions for this but I have some trouble with the computeTextureCoordinates part. Since quantities are now computed on faces, we need to make loops around vertices on the dual mesh. The problem is that in the original Stripe Patterns paper the "frequency adjustment" step (p.7) is only defined for triangles and I'm not sure if it can be generalized to n-gons (maybe @keenancrane has an idea about this?)

[nmwsharp]

Hi @DavidJourdan thanks again for this! I went ahead and merged this since the code looks pretty polished :) Hopefully it can be useful to others.

However, I tried running it in a demo app to take some screenshots for docs, and couldn't get it to work. Can you tell what I'm doing wrong? Here's the the important snippet of the code I tried (it's set up to run the latest commits of Polyscope and geometry-central).

It's possible that there are some Polyscope problems going on visualizing the corner data & direction field, I made some recent changes to that stuff in Polyscope. But that wouldn't effect the isolines curves, and they seem off too.

void myCallback() {

  ImGui::SliderFloat("param", &freq, 0., 100.);

  if (ImGui::Button("do work")) {

    // Get vertex tangent spaces
    geometry->requireVertexTangentBasis();
    VertexData<Vector3> vBasisX(*mesh);
    VertexData<Vector3> vBasisY(*mesh);
    for (Vertex v : mesh->vertices()) {
      vBasisX[v] = geometry->vertexTangentBasis[v][0];
      vBasisY[v] = geometry->vertexTangentBasis[v][1];
    }

    VertexData<Vector2> guideField =
        geometrycentral::surface::computeSmoothestVertexDirectionField(
            *geometry);
    psMesh->addVertexTangentVectorQuantity("guide field", guideField, vBasisX,
                                           vBasisY, 2);

    VertexData<double> frequencies(*mesh, freq);

    CornerData<double> periodicFunc;
    FaceData<int> zeroIndices;
    FaceData<int> branchIndices;

    std::tie(periodicFunc, zeroIndices, branchIndices) =
        computeStripePattern(*geometry, frequencies, guideField);

    for (Corner c : mesh->corners()) {
      periodicFunc[c] = std::fmod(periodicFunc[c], 2 * M_PI);
    }

    psMesh->addCornerScalarQuantity("periodic", periodicFunc);

    std::vector<Vector3> isolineVerts;
    std::vector<std::array<size_t, 2>> isolineEdges;
    std::tie(isolineVerts, isolineEdges) = extractPolylinesFromStripePattern(
        *geometry, periodicFunc, zeroIndices, branchIndices, guideField, false);

    polyscope::registerCurveNetwork("isolines", isolineVerts, isolineEdges);
  }
}

int main(int argc, char **argv) {

  // Configure the argument parser
  args::ArgumentParser parser("geometry-central & Polyscope example project");
  args::Positional<std::string> inputFilename(parser, "mesh", "A mesh file.");

  // Parse args
  try {
    parser.ParseCLI(argc, argv);
  } catch (args::Help &h) {
    std::cout << parser;
    return 0;
  } catch (args::ParseError &e) {
    std::cerr << e.what() << std::endl;
    std::cerr << parser;
    return 1;
  }

  // Make sure a mesh name was given
  if (!inputFilename) {
    std::cerr << "Please specify a mesh file as argument" << std::endl;
    return EXIT_FAILURE;
  }

  // Initialize polyscope
  polyscope::init();

  // Set the callback function
  polyscope::state::userCallback = myCallback;

  // Load mesh
  std::tie(mesh, geometry) = readManifoldSurfaceMesh(args::get(inputFilename));

  // Register the mesh with polyscope
  psMesh = polyscope::registerSurfaceMesh(
      polyscope::guessNiceNameFromPath(args::get(inputFilename)),
      geometry->inputVertexPositions, mesh->getFaceVertexList());

  psMesh->setAllPermutations(polyscopePermutations(*mesh));

  // Give control to the polyscope gui
  polyscope::show();

  return EXIT_SUCCESS;
}

[DavidJourdan]

Hi, thanks for merging this!

From what I can see, the issue seems to be with the fmod rounding (and also the guideField which should be a 2-RoSy). This modified version should work:

void myCallback() {

  ImGui::SliderFloat("param", &freq, 0., 100.);

  if (ImGui::Button("do work")) {

    // Get vertex tangent spaces
    geometry->requireVertexTangentBasis();
    VertexData<Vector3> vBasisX(*mesh);
    VertexData<Vector3> vBasisY(*mesh);
    for (Vertex v : mesh->vertices()) {
      vBasisX[v] = geometry->vertexTangentBasis[v][0];
      vBasisY[v] = geometry->vertexTangentBasis[v][1];
    }

    VertexData<Vector2> guideField =
        geometrycentral::surface::computeSmoothestVertexDirectionField(
            *geometry, 2);
    psMesh->addVertexTangentVectorQuantity("guide field", guideField, vBasisX,
                                           vBasisY, 2);

    VertexData<double> frequencies(*mesh, freq);

    CornerData<double> periodicFunc;
    FaceData<int> zeroIndices;
    FaceData<int> branchIndices;

    std::tie(periodicFunc, zeroIndices, branchIndices) =
        computeStripePattern(*geometry, frequencies, guideField);

    // for (Corner c : mesh->corners()) {
    //   periodicFunc[c] = std::fmod(periodicFunc[c], 2 * M_PI);
    // }

    psMesh->addCornerScalarQuantity("periodic", periodicFunc);

    std::vector<Vector3> isolineVerts;
    std::vector<std::array<size_t, 2>> isolineEdges;
    std::tie(isolineVerts, isolineEdges) = extractPolylinesFromStripePattern(
        *geometry, periodicFunc, zeroIndices, branchIndices, guideField, false);

    polyscope::registerCurveNetwork("isolines", isolineVerts, isolineEdges);
  }
}

int main(int argc, char **argv) {

  // Configure the argument parser
  args::ArgumentParser parser("geometry-central & Polyscope example project");
  args::Positional<std::string> inputFilename(parser, "mesh", "A mesh file.");

  // Parse args
  try {
    parser.ParseCLI(argc, argv);
  } catch (args::Help &h) {
    std::cout << parser;
    return 0;
  } catch (args::ParseError &e) {
    std::cerr << e.what() << std::endl;
    std::cerr << parser;
    return 1;
  }

  // Make sure a mesh name was given
  if (!inputFilename) {
    std::cerr << "Please specify a mesh file as argument" << std::endl;
    return EXIT_FAILURE;
  }

  // Initialize polyscope
  polyscope::init();

  // Set the callback function
  polyscope::state::userCallback = myCallback;

  // Load mesh
  std::tie(mesh, geometry) = readManifoldSurfaceMesh(args::get(inputFilename));

  // Register the mesh with polyscope
  psMesh = polyscope::registerSurfaceMesh(
      polyscope::guessNiceNameFromPath(args::get(inputFilename)),
      geometry->inputVertexPositions, mesh->getFaceVertexList());

  psMesh->setAllPermutations(polyscopePermutations(*mesh));

  // Give control to the polyscope gui
  polyscope::show();

  return EXIT_SUCCESS;
}

[nmwsharp]

Oh yeesh, not using the 2 RoSy field was a bad mistake! Thanks for pointing that out, I can confirm that it's producing much more sane outputs with that fix :) I'll pull a quick docs page together.

[nmwsharp]

Now live with docs here! https://geometry-central.net/surface/algorithms/stripes/ Feel free to submit changes if anything is needed.