stack overflow in insertEdges()

[msokalski]

For few not so trivial cases (actually pretty extreme ones), I'm getting stack overflow exception from insertEdge() call. There are many many CDT::Triangulation<double,CDT::LocatorKDTree<double,32,32,32>>::triangulatePseudopolygon calls on the stack causing this. That would be great if refactoring this recursion into a pseudo-recursion loop would be possible. I'm on Win10, VS22. Problem should be reproducible with something like:

#include <random>
#include "CDT.h"

#define POINTS 1000000
#define EDGES 10000
std::random_device rd{};
std::mt19937_64 gen{rd()};
std::gamma_distribution<double> d(0.1,2.0);
std::vector<CDT::V2d<double>> verts;
for (int i = 0; i < POINTS; i++)
{
  CDT::V2d<double> v;
  v.x = d(gen) - 5.0;
  v.y = d(gen) - 5.0;
  verts.push_back(v);
}

std::vector<CDT::Edge> edges;
for (int c = 0; c < EDGES; c++)
{
  int from = rand()%POINTS;
  int to = (from + 1 + rand() % (POINTS - 1)) % POINTS;
  CDT::Edge e{ (size_t)from, (size_t)to };
  edges.push_back(e);
}

CDT::RemoveDuplicatesAndRemapEdges(verts,edges);

CDT::Triangulation<double> cdt;
cdt.insertVertices(verts);
cdt.insertEdges(edges);  // here it occurs
cdt.eraseSuperTriangle();

EDIT: Enlarging stack reserve size in linker (from a default of 1MB) to 10MB resolves the issue (for 1M verts and 10K edges).

[artem-ogre]

Could it be caused by duplicated points or overlapping edges? CDT does not support neither.

From README:

Pre-conditions:

• No duplicated points (use provided functions for removing duplicate points and re-mapping edges)
• No two constraint edges intersect each other (overlapping boundaries are allowed)

If thats not the case, could you provide a concrete input configuration that can be used to deterministically reproduce the crash?

I’m 95% sure it’s duplicates or overlaps though. Issues like this get raised periodically up to a point that I should probably emphasize it more in documentation and provide issue creation guidelines.

[artem-ogre]

For issues caused by wrong input see #69 #24 #2 or https://github.com/artem-ogre/CDT/issues?q=label%3A%22wrong+input%22+is%3Aclosed

[msokalski]

Hi @artem-ogre , No it is not the input problem, in the given example it is pre-filtered using RemoveDuplicatesAndRemapEdges. Also this is not CDT bug, but only refactoring suggestion. CDT is easily able to solve this particular input but it requires changing linking parameters.

Currently insertEdges() makes use of recursive pseudo-polygon algorithm. It is elegant and is fast but it is fragile if too small stack size is given at the compile/link time for executable.

Of course, one could increase executable stack size if scale of the problem can be predicted at the compile time, but usually it is not possible.

One way to resolve this uncomfortable situation is replacing recursive calls to 'triangulatePseudopolygon' with a loop emulating call stack (arguments, local variables) on the heap rather on stack.

I hope this time I'm more clear :)

[artem-ogre]

Yes, I understand. I could re-write recursion into stack-based approach. However it is hard for me to imagine the kind of a pseudo-poly that is so big that it causes the stack overflow.

What about crossing edges? I don't see any checks or fixes for crossing edges. Could it be that this is what causing the pseudo-polygon that huge.

[msokalski]

I've seen CDT was able to 'survive' all 3 cases: crossing edges, edge-edge overlapping and edge overlapping multiple vertices. This is record breaking! I think in this particular example, pseudo polygons are very asymmetric, so splitting them is far from halving what makes stack size requirement more like N than log(N).

[artem-ogre]

Edge-edge overlapping and edge overlapping multiple vertices are supported. And there is quite some internal complexity involved to correctly handle such corner cases :) There is also a branch that implements automatically resolving crossing edges intersection points (not 100% numerically robust as intersection points are subject to rounding errors but works in most practical cases).

Surviving crossing edges is not a big deal, the problem is that there are no guarantees that the output is correct.

[msokalski]

Ok, I'm going to add edge-crossing check in order to see if that would reduce stack size requirement. I'll let you know after the weekend.

[msokalski]

Hi, good news: resolving crossing edges greatly reduces stack size needed for constraining. Indeed, max number of polygons edges is reduced 10x or so. I've also found out that vertex indices provided to edges vector for insertEdges().must be decremented by 3, after calling RemoveDuplicatesAndRemapEdges(). Is this intentional? I believe it has something to do with a super-triangle, but I feel a bit confused. In example if I want to add edge over indices: {4,1} I need to provide {1,(size_t)-2} a negative index!

[artem-ogre]

Could you provide a simple example of the indexing problem? There should be absolutely no index mangling required from the user.

Indeed super-triangle will offset edge indices by 3, but this should be taken care of internally, when calling insertEdges. You can see in the following code that 3 will be added to vertex indices of edge (m_nTargetVerts is 3):

template <typename T, typename TNearPointLocator>
template <
    typename TEdgeIter,
    typename TGetEdgeVertexStart,
    typename TGetEdgeVertexEnd>
void Triangulation<T, TNearPointLocator>::insertEdges(
    TEdgeIter first,
    const TEdgeIter last,
    TGetEdgeVertexStart getStart,
    TGetEdgeVertexEnd getEnd)
{
    for(; first != last; ++first)
    {
        // +3 to account for super-triangle vertices
        const Edge edge(
            VertInd(getStart(*first) + m_nTargetVerts),
            VertInd(getEnd(*first) + m_nTargetVerts));
        insertEdge(edge, edge);
    }
    eraseDummies();
}

After removing the super-triangle with e.g., eraseSuperTriangle, edges will be mapped back. The following code will be called:

inline Edge RemapNoSuperTriangle(const Edge& e)
{
    return Edge(e.v1() - 3, e.v2() - 3);
}

[msokalski]

Without enabling fix_indices, app just crashes.

    {
        bool fix_indices = false;

        struct MyPoint
        {
            MyPoint(double x, double y) : x(x), y(y) {}
            double x, y;
        };

        std::vector<MyPoint> cloud;
        cloud.push_back(MyPoint(-10, -10)); //    2--___3
        cloud.push_back(MyPoint(+10, -10)); //    |     |
        cloud.push_back(MyPoint(-10, +10)); //    |      |
        cloud.push_back(MyPoint(+9, +9));   //    0------1  

        struct MyEdge
        {
            MyEdge(int a, int b) : a(a), b(b) {}
            int a, b;
        };

        std::vector<MyEdge> constr;
        constr.push_back(MyEdge(1,2)); // force longer diagonal

        std::vector<CDT::V2d<double>> nodups;
        for (int i = 0; i < cloud.size(); i++)
        {
            CDT::V2d<double> v;
            v.x = cloud[i].x;
            v.y = cloud[i].y;
            nodups.push_back(v);
        }

        std::vector<CDT::Edge> edges;
        for (int c = 0; c < constr.size(); c++)
        {
            CDT::Edge e
            {
                (size_t)constr[c].a,
                (size_t)constr[c].b
            };
            edges.push_back(e);
        }

        CDT::DuplicatesInfo dups_nfo = CDT::RemoveDuplicatesAndRemapEdges(nodups, edges);

        // magic part
        if (fix_indices)
        {
            int account_super_triangle = 3;
            for (int c = 0; c < edges.size(); c++)
            {
                CDT::Edge e = edges[c];
                CDT::Edge f{ e.v1() - account_super_triangle,e.v2() - account_super_triangle };
                edges[c] = f;
            }
        }

        CDT::Triangulation<double> cdt;
        cdt.insertVertices(nodups);
        cdt.eraseSuperTriangle();

        cdt.insertEdges(edges);

        printf("triangles: %d\n", cdt.triangles.size());
        for (int i = 0; i < cdt.triangles.size(); i++)
        {
            printf("%d %d %d\n",
                cdt.triangles[i].vertices[0],
                cdt.triangles[i].vertices[1],
                cdt.triangles[i].vertices[2]);
        }

        exit(0);
    }

[artem-ogre]

You should not call insertEdges after eraseSuperTriangle: after eraseSuperTriangle super-triangle vertices are erased, but +3 adjustment will still be applied inside insertEdges. Generally, one should think of eraseXXX methods as a final 'bake-the-final-data' step.

[artem-ogre]

This is not a great API design and especially not very newcomer-friendly.

[msokalski]

Oh great! I didn't realize calling eraseSuperTriange modifies vertex indices (I thought it affects triangles only). Thanks for your help!

[msokalski]

FYI, trying to add an edge over many colinear vertices which are already edge connected (like in the test case for https://github.com/artem-ogre/CDT/issues/84) can easily overflow stack in the following recursion: https://github.com/artem-ogre/CDT/blob/21cdbb249e58e5628a61ac1a0a1753ee0ec5d437/CDT/include/CDT.hpp#L469

I can see there are many commits made since I pulled, so this may be outdated :)

[artem-ogre]

There is always an input big enough to overflow the stack :) The question is: does this happen in real world with the real data or does this happen with specifically designed synthetic cases only? Is it causing real pain or is it mostly of academic interest?

If it is the latter, I should (probably) fix this anyway, but it will be very low on priority and may not happen in the foreseeable future.

[msokalski]

I can barely imagine such real world input but who knows :) If you feel this is hard to reimplement in recursion free way or it would made code significantly harder to maintain, I'd suggest adding depth counter to CDT::Triangulation and a user settable depth limit, so you could throw a regular C++ exception instead of dying with SO. I use CDT for validating corner cases in my own lib, so yes, you can keep it as a bottom priority thing, I'm not real world user :) Anyway, thank you so much for making such beautiful code, open source!

[artem-ogre]

You should not call insertEdges after eraseSuperTriangle: after eraseSuperTriangle super-triangle vertices are erased, but +3 adjustment will still be applied inside insertEdges. Generally, one should think of eraseXXX methods as a final 'bake-the-final-data' step.

This is not a great API design and especially not very newcomer-friendly.

@msokalski I have changed the code so that it would detect if erase... was called (see new Triangulation::isFinalized method) and would throw an exception when trying to add new vertices or edges. This should make the API more error-proof.

[artem-ogre]

Hello, @msokalski, I converted code for inserting/conforming-to edges from recursion to iteration (see the branch). Could you please check if it addresses the problems you had?

[msokalski]

Yes, sure. Ill check it on Monday, I'm away for weekend. Thanks.

[artem-ogre]

great, thanks! have a nice weekend :)

[msokalski]

Hey @artem-ogre , that was a long weekend :p Actually I was quite busy with my own stuff, so please accept my apologies for such delay!

I've confirmed, no more stack overflows occurs in insertEdges(), I didn't test conforming edges as it is currently outside my scope, but probably it is also fine (I've reviewed sources while bench marking). I've noticed a tiny performance drop, about 5% when compared with stack based implementation but I still prefer stability over speed so thanks for changing it!

[artem-ogre]

It is totally fine, it is your free time and I really appreciate it. Thank you!

[artem-ogre]

Solved by #102