Missing element in triangulation of point cloud due to small super-triangle

[AndreFecteau]

Recently I have been evaluating constrained Delaunay triangulation libraries and have stumbled onto this gem. During my testing of this tool I have encountered one triangulation where I believe is not the desired output but due to having the super-triangle being too small for the point cloud. Here is the point cloud that causes the issue also attached you can find the visualization of the triangulation including the super-triangle and missing edge in red.

PointCloudTriangulation.pdf

VERTICES
0.15147567518991,-5.144230291488,0
0.10442,-5.06098,0
0.1228735248885,-5.12897066233,0
0.1322969035965,-5.141286788425,0
0.115882234089,-5.115648852375,0
0.13262891777369,-5.119598039304,0
0.11311813200326,-5.138343285364,0
0.111139345548,-5.09040217055,0
0.121648346645,-5.09478239621,0
0.1152696449675,-5.098554719315,0
0.1303134549875,-5.106097900345,0
0.10889094329024,-5.1023270424231,0
0.142052296482,-5.131914165395,0
0.154751587595,-5.127544895745,0
0.1484019420385,-5.129729530575,0
0.1389785633305,-5.117413404475,0
0.1563895437975,-5.119202197875,0
0.150039898241,-5.1213868327,0
0.140515429906,-5.12466378494,0
0.145328208887,-5.11522876965,0
0.1516778544435,-5.113044134825,0
0.1580275,-5.1108595,0
0.107779672774,-5.075691085275,0
0.112085873903,-5.07104664934,0
0.1514750401785,-5.061655309135,0
0.13440575,-5.08723775,0
0.110784,-5.063616,0
0.11338774780615,-5.078477298681,0
0.1682004187975,-5.131013072875,0
0.18164925,-5.13448125,0
0.1921660803575,-5.05969461827,0
0.1750967901785,-5.085277059135,0
0.221046245536,-5.045923052405,0
0.31820357142857,-4.9298217230895,0
0.2499264107145,-5.032151486545,0
0.284064991072,-4.980986604815,0

[artem-ogre]

Hello, @AndreFecteau

Thank you very much for the issue and a detailed description with an example! :)

If any vertex belongs to the super-triangle during edge-flip test: algorithm should treat this as a special case and handle it appropriately. So the size of the super-triangle should not matter.

I will check if there is a bug in that special-case handling.

[artem-ogre]

The bug was there indeed, thank you again :)

Please check this PR: #38
The fix is in the branch bugfix/37-fix-flip-needed-test.
If the fix works for your datasets, I will add a test file, do some more testing and merge it.

[AndreFecteau]

Thanks for the quick fix, I will try this in the next few days.

[AndreFecteau]

I ran it through all the benchmarks and the fix seemed to have worked. Thanks a lot!!

[artem-ogre]

@AndreFecteau May I ask what did you compare CDT with and how did it compare? This is very interesting :)

[AndreFecteau]

@artem-ogre These numbers are only meant to be a rough measurement. Note that not all these libraries are reliable or support general constrained Delaunay triangulation.


std::vector<std::pair<double, double>> vertices;
for(size_t i = 0; i < 1000; ++i)
{
    for(size_t j = 0; j < 1000; ++j)
    {
        vertices.emplace_back(i*1e-5, j*1e-5);
    }
}


10x10

Elapsed time poly2Tri: 0.00012235 s
Elapsed time CDT: 0.000653747 s
Elapsed time CGAL: 0.000221772 s
Elapsed time Delaunator: 6.9119e-05 s
Elapsed time Triangle: 0.000959387 s

20x20

Elapsed time poly2Tri: 0.000755585 s
Elapsed time CDT: 0.00217176 s
Elapsed time CGAL: 0.000819186 s
Elapsed time Delaunator: 0.000294325 s
Elapsed time Triangle: 0.00246849 s


100x100

Elapsed time poly2Tri: 0.0116399 s
Elapsed time CDT: 0.113256 s
Elapsed time CGAL: 0.0159775 s
Elapsed time Delaunator: 0.00511762 s
Elapsed time Triangle: 0.0100751 s

100x1000

Elapsed time poly2Tri: 0.106778 s
Elapsed time CDT: 8.49504 s
Elapsed time CGAL: 0.158438 s
Elapsed time Delaunator: 0.0808353 s
Elapsed time Triangle: 0.127281 s

1000x1000

Elapsed time poly2Tri: 1.1229 s
Elapsed time CDT: 84.6357 s
Elapsed time CGAL: 1.6105 s
Elapsed time Delaunator: 0.964066 s
Elapsed time Triangle: 1.69101 s




srand (time(NULL));
std::vector<std::pair<double, double>> vertices;
for(size_t i = 0; i < 10; ++i)
{
    for(size_t j = 0; j < 10; ++j)
    {
        bool inserted = false;
        while(!inserted)
        {
            size_t x = (rand() % 10000);
            size_t y = (rand() % 10000);
            auto it = std::find_if(vertices.begin(), vertices.end(), [&](const auto& xy){return x-xy.first < 1e-5 && y-xy.second < 1e-5;});
            if(it == vertices.end())
            {
                vertices.emplace_back(x*1e-5, y*1e-5);
                inserted = true;
            }
        }
    }
}

10x10

Elapsed time poly2Tri: 0.000154201 s
Elapsed time CDT: 0.000439808 s
Elapsed time CGAL: 9.577e-05 s
Elapsed time Delaunator: 6.046e-05 s
Elapsed time Triangle: 0.000652088 s



100x100

Elapsed time poly2Tri: 0.0123536 s
Elapsed time CDT: 0.0613802 s
Elapsed time CGAL: 0.00735314 s
Elapsed time Delaunator: 0.00525829 s
Elapsed time Triangle: 0.0123408 s


100x1000

Elapsed time poly2Tri: 0.168798 s
Elapsed time CDT: 1.06543 s
Elapsed time CGAL: 0.0847132 s
Elapsed time Delaunator: 0.0803343 s
Elapsed time Triangle: 0.191198 s


500x500

Elapsed time poly2Tri: 0.448719 s
Elapsed time CDT: 3.0024 s
Elapsed time CGAL: 0.213624 s
Elapsed time Delaunator: 0.238802 s
Elapsed time Triangle: 0.560273 s

[artem-ogre]

@AndreFecteau Is CDT using boost with CDT_USE_BOOST in these benchmarks?

[artem-ogre]

And if yes, does it use nearestVertexRtree?

[AndreFecteau]

Yes CDT_USE_BOOST is enabled, I attached a flame graph of the benchmark 500x500. While I don't have experience with nearest neighbour search in point tree, I've always found that Boost::RTree access iterators are quite expensive for simple shapes. Have you tried using a Oct Tree or kd-Tree for this section of the code?

[artem-ogre]

Really appreciate the input, thank you! I will look into it as soon as I get some free time. I will open a new issue for this.

[artem-ogre]

Yes CDT_USE_BOOST is enabled, I attached a flame graph of the benchmark 500x500. While I don't have experience with nearest neighbour search in point tree, I've always found that Boost::RTree access iterators are quite expensive for simple shapes. Have you tried using a Oct Tree or kd-Tree for this section of the code?

@AndreFecteau
I got some time to do some profiling in "release with debug info" configuration. Even though I get similar timings to yours, flame graph looks very different and r-tree performance is not a bottleneck (see more info below, original interactive flame graph:flamegraph.zip). Do you have any thoughts on that? What was your benchmarked configuration? If you answer could you please answer in #40?

Flame Graph

Benchmarked code

#include "CDT.h"

typedef double CoordType;
typedef CDT::Triangulation<CoordType> Triangulation;
typedef CDT::V2d<CoordType> V2d;
typedef CDT::Vertex<CoordType> Vertex;
typedef CDT::Triangle Triangle;
typedef CDT::Box2d<CoordType> Box2d;
typedef CDT::Edge Edge;

#include <chrono>
#include <iostream>

int main(int argc, char* argv[])
{
    const size_t N = 1000;
    std::vector<V2d> vertices;
    vertices.reserve(N * N);
    for(size_t i = 0; i < N; ++i)
    {
        for(size_t j = 0; j < N; ++j)
        {
            vertices.push_back(V2d::make(i * 1e-5, j * 1e-5));
        }
    }

    using std::chrono::duration;
    using std::chrono::duration_cast;
    using std::chrono::high_resolution_clock;
    using std::chrono::milliseconds;
    auto t1 = high_resolution_clock::now();

    Triangulation triangulation(CDT::FindingClosestPoint::BoostRTree);
    triangulation.insertVertices(vertices);

    auto t2 = high_resolution_clock::now();
    /* Getting number of milliseconds as an integer. */
    auto ms_int = duration_cast<milliseconds>(t2 - t1);
    std::cout << ms_int.count() << "ms\n";

    return 0;
}