How to update triangulations efficiently? about cdt

[artem-ogre]

Could you elaborate more on what do you aim to achieve by splitting edges? Why do you need to split them and based on what criteria?

[pan3rock]

I'm trying to follow the work SAGRPF. The algorithm mainly described some criteria to determine which edges to split in the next iteration. I transform the algorithm to my work but face the problem the cost of constructing triangulation is too high. Due to the fact that only a few new points are added at each iteration, is there some way to reduce the time?

I run three examples, and the cost of the constructing is below:

the "mesh" timing is this part

    CDT::Triangulation<double> cdt;
    cdt.insertVertices(nodes_);
    cdt.eraseSuperTriangle();
    auto elements = cdt.triangles;
    auto edges = CDT::extractEdgesFromTriangles(elements);

Now I'm trapped in the situation that the cost of evaluating the target function is much less than that of the algorithm, especially the construction of triangulation, so scanning directly in finer grid may be better than through the algorithm.

[artem-ogre]

You could try making Triangulation::splitFixedEdgeAt method public and calling it directly to avoid re-triangulating everything:

CDT/CDT/include/Triangulation.h

Line 684 in 58f34da

VertInd splitFixedEdgeAt(

[artem-ogre]

Also you might want to expose Triangulation::edgeTriangles method to get the triangles sharing an edge.

[pan3rock]

I tried to follow your advice but failed. I create a simple demo as below

#include "timer.hpp"
#include <CDT.h>
#include <Eigen/Dense>
#include <fmt/format.h>
#include <highfive/H5Easy.hpp>
#include <iostream>
#include <vector>

using namespace Eigen;
using V2d = CDT::V2d<double>;

Eigen::ArrayXd compute_edge_length(const std::vector<V2d> &nodes,
                                   const CDT::EdgeVec &edges) {
  std::vector<double> le;
  for (auto it = edges.begin(); it != edges.end(); ++it) {
    auto v1 = nodes[it->v1()];
    auto v2 = nodes[it->v2()];
    double d = CDT::distance(v1, v2);
    le.push_back(d);
  }
  ArrayXd ret = Map<ArrayXd, Unaligned>(le.data(), le.size());
  return ret;
}

int main() {
  std::vector<V2d> nodes_append;
  std::vector<double> x_init{0, 1, 1, 0}, y_init{0, 0, 1, 1};
  for (int i = 0; i < 4; ++i) {
    nodes_append.push_back(V2d::make(x_init[i], y_init[i]));
  }

  std::vector<V2d> nodes;
  const int maxiter = 100;
  Timer::begin("mesh");
  CDT::Triangulation<double> cdt;
  cdt.insertVertices(nodes_append);
  cdt.eraseSuperTriangle();
  Timer::end("mesh");

  for (int i = 0; i < maxiter; ++i) {
    for (auto it = nodes_append.begin(); it != nodes_append.end(); ++it) {
      nodes.push_back(*it);
    }

    Timer::begin("mesh");
    auto elements = cdt.triangles;
    auto edges_us = CDT::extractEdgesFromTriangles(elements);
    Timer::end("mesh");

    CDT::EdgeVec edges(edges_us.begin(), edges_us.end());

    // for simplicity, choose the longest edge to split
    ArrayXd edge_length = compute_edge_length(nodes, edges);
    int imax;
    edge_length.maxCoeff(&imax);

    auto edge_lmax = edges[imax];
    auto n1 = nodes[edge_lmax.v1()];
    auto n2 = nodes[edge_lmax.v2()];
    nodes_append.clear();
    auto na = V2d::make((n1.x + n2.x) / 2.0, (n1.y + n2.y) / 2.0);
    nodes_append.push_back(na);

    Timer::begin("mesh");
    auto inb = cdt.edgeTriangles(edge_lmax.v1(), edge_lmax.v2());
    auto isplit = cdt.splitFixedEdgeAt(edge_lmax, na, inb.first, inb.second);
    Timer::end("mesh");
    fmt::print("{:10d}{:10d}\n", i, isplit);
  }

  const int num_node = nodes.size();
  ArrayXd x(num_node), y(num_node);
  for (int i = 0; i < num_node; ++i) {
    x(i) = nodes[i].x;
    y(i) = nodes[i].y;
  }

  auto elements = cdt.triangles;
  ArrayXXi tri(elements.size(), 3);
  for (size_t i = 0; i < elements.size(); ++i) {
    tri(i, 0) = elements[i].vertices[0];
    tri(i, 1) = elements[i].vertices[1];
    tri(i, 2) = elements[i].vertices[2];
  }

  H5Easy::File fout("demo2.h5", H5Easy::File::Overwrite);
  H5Easy::dump(fout, "x", x);
  H5Easy::dump(fout, "y", y);
  H5Easy::dump(fout, "tri", tri);

  std::cout << Timer::summery() << std::endl;

  return 0;
}

the output errors

It seems m_vertTris is empty.

Though I can write it by myself through unordered_map, there is a problem: when the edge is on the boundary, in other words, there is only one triangle adjacent to the edge, what to return?

Maybe it's too difficult for me to account for speeding it up.

[artem-ogre]

You can avoid boundary edges by not calling eraseSuperTriangle (do it at the very end instead) and not touching the edges where both vertices have index less than 3 (first three vertices are super-triangle vertices).

[pan3rock]

Does it mean that I need to subtract 3 from index when using vertices? e.g.

auto v1 = edge.v1()
auto v2 = edge.v2()
double x = nodes[v1 - 3].x;

except super ones (0, 1, 2) the equation is satisifed:
Index of vertices(with eraseSuperTriangle) = Index of vertices(without eraseSuperTriangle) - 3
is that right?

I've completed the modification, and it seems good at the simple example (x4 faster). I'll test it in my project.

[artem-ogre]

Yes, you're right. If you don't call eraseSuperTriangle you need to add that offset of 3.

But I also meant that you should never split super-triangle edges:

if(edge.v2() < 3) // edge vertices are sorted (v1 <= v2), so only a single comparison is needed 
{
    // skip the edge, don't calculate it's length
}

[pan3rock]

When the number of splitted edges for each iteration are greater than 1, edgeTriangles failed. It seems that splitFixedEdgeAt changed the original triangles, so the second splitted edges can not be found.

codes below:

for (int i = 0; i < maxiter; ++i) {
  ...
  edges_split = ... // greater than 1
  for (auto it = edges_split_.begin(); it != edges_split_.end(); ++it) {
        auto n1 = nodes_[it->v1()];
        auto n2 = nodes_[it->v2()];
        double x = (n1.x + n2.x) / 2.0;
        double y = (n1.y + n2.y) / 2.0;
        auto na = V2d::make(x, y);
        auto inb = cdt_.edgeTriangles(it->v1(), it->v2());
        auto m = cdt_.splitFixedEdgeAt(*it, na, inb.first, inb.second);
  }
  ...
}

The initial triangles:

add the point from the first edge_split (0, 0.5), the original triangles have been changed.

So the second point (0.5, 0.5) can't be added.

[pan3rock]

Thanks again. I found it's not gonna work in this way. Mesh construction has always accounted for the bulk of time in mesh-dependent algorithms, e.g. Sambridge's neighborhood algorithm (voronoi). I'll try other ideas.