Degenerated triangulations, how to extract polygons?

[msokalski]

Is this possible, without re-inventing the wheel, to detect and extract info on degenerated triangulations?
Simply grouping triangles into 'polygons' would be sufficient in my case, although extracting N-gon mesh (with adjacencies, etc) is the nicest way to go.

    template <typename T>
    std::vector<CDT::V2d<T>> genData(const std::size_t size)
    {
        const T x = 0x5af2efc1.p-30;
        const T y = 0x348268e0.p-30;
        const T r = 0x6904d1c1.p-30;

        const T dx = 2 * x;
        const T dy = 2 * (r + y);

        const CDT::V2d<T> p[4] =
        {
            CDT::V2d<T>::make(0,0),
            CDT::V2d<T>::make(0,2 * y),
            CDT::V2d<T>::make(x,r + y),
            CDT::V2d<T>::make(x,r + 3 * y),
        };

        size_t rows = (size_t)ceil(sqrt(size / 7.0));
        size_t cols = (size + 2 * rows + 1) / (4 * rows);

        std::vector<CDT::V2d<T>> vv;
        vv.reserve(4 * cols * rows);

        for (int row = 0; row < rows; row++)
            for (int col = 0; col < cols; col++)
                for (int i = 0; i < 4; i++)
                    vv.push_back(CDT::V2d<T>::make(
                        p[i].x + col * dx, 
                        p[i].y + row * dy));

        return vv;
    }

[artem-ogre]

Not sure I understand what is shown in the image and what is it you would like to achieve?

[msokalski]

Ok, I'm sorry being not clear. The image presents 2 valid triangulations of exactly same input, red and green. Yellow edges are ones that are identical in both triangulations. What I would like to achieve is a way to extract 'unified polygonal mesh' for degenerated parts of triangulations (triangles adjacent to each other with vertices sharing same circumcenter).

[artem-ogre]

I am red-green colorblind, so this is really hard to read for me :) (usually cyan, yellow and magenta are safe picks in such cases)

[artem-ogre]

Would you like to classify each triangle by the hexagon it belongs to?
If you can calculate hexagons that vertices belong to, you can then for each triangle choose the one and only hexagon that all 3 triangle vertices share.

[msokalski]

Ok, here 2 triangulations are: yellow and cyan. White color indicates common solution.

I hope it's ok now.

[msokalski]

Hexagons is just simple example, I'd like to be able to extract such polygons for any input given where I don't really know where such degenerations would occur.

[artem-ogre]

This case is not really a degeneration: this happens because of points lie exactly on the circumcircles of some adjacent triangles so there is more than one way to connect points into a Delaunay triangulation.

Imagine a square: (0,0), (1,0), (1,1), (0, 1). In this case both ways to draw a diagonal will produce a triangulation satisfying Delaunay property: "A circle circumscribing any Delaunay triangle does not contain any other input points in its interior."

So for you to detect such spots, you could go over each triangle and for the 'opposite' vertex of each edge do the robust in-circle test and when it is exactly 0 you have found a case where edge flip still maintains Delaunay property.

[msokalski]

Is this possible to make such tests using current CDT api or should I route my code directly to predicates library?

[artem-ogre]

Did not quite understand all the details, but looks okay to me.
I've done some mostly cosmetic fixes and modernized the snippet. Also replaced list with vector. I did not test it, so it very well might be broken.

typedef std::vector<TriInd> Poly;
template <typename T, typename L = LocatorKDTree<T> >
std::vector<Poly> Polygonize(const Triangulation<T, L>& cdt)
{
    typedef TriInd PolyInd;
    // fer every tri, here we store index of the poly the tri belongs to
    auto map = std::vector<PolyInd>(cdt.triangles.size());
    // vector of Polys, each as a list of tri indices belonging to the Poly
    auto polys = std::vector<Poly>();

    for(TriInd iT = 0; iT < cdt.triangles.size(); iT++)
    {
        const auto& tri = cdt.triangles[iT];
        auto merged = false;

        // compare i'th tri with its adjacent tris
        for(const auto adj : tri.neighbors)
        {
            // but only if adjacent is processed already
            if(adj > iT)
                continue;
            // locate reflex vert in adj triangle
            const auto& vr =
                cdt.vertices[opposedVertex(cdt.triangles[adj], iT)];
            const auto& v1 = cdt.vertices[tri.vertices[0]];
            const auto& v2 = cdt.vertices[tri.vertices[1]];
            const auto& v3 = cdt.vertices[tri.vertices[2]];
            using predicates::adaptive::incircle;
            if(!incircle(vr.x, vr.y, v1.x, v1.y, v2.x, v2.y, v3.x, v3.y))
            {
                if(!merged)
                {
                    // append tri to already existing poly
                    merged = true;
                    const auto append_to = map[adj];
                    map[iT] = append_to;
                    polys[append_to].push_back(iT);
                }
                else
                {
                    const auto merge_to = map[iT];
                    const auto merge_from = map[adj];

                    if(merge_to == merge_from)
                        continue;

                    // funny case, tri is a bridge between 2 polys merge'em all
                    // together
                    for(const auto i : polys[merge_from])
                    {
                        map[i] = merge_to;            // remap
                        polys[merge_to].push_back(i); // merge
                    }
                    if(merge_from != polys.size() - 1)
                    {
                        // replace merge_from poly with last poly in polys
                        polys[merge_from] = polys.back();
                        for(const auto i : polys[merge_from])
                        {
                            map[i] = merge_from; // remap
                        }
                    }
                    polys.pop_back();
                }
            }
        }

        if(!merged)
        {
            // at the moment, just alone tri
            // make a new poly for it
            map[iT] = polys.size();
            polys.push_back({iT});
        }
    }

    return polys;
}

[msokalski]

Hey @artem-ogre
Looks way more like a part of CDT now :)
The most I like is refactored for loops!
Thank you.

[msokalski]

Hi @artem-ogre,
I'm adding post processing to this routine, allowing user to easily extract every polygon's boundary as a sequence of vertices.
To achieve this, I need to "rotate" triangles, more specifically ensure that specified vertex of a triangle appears on Triangle::vertices[0], also I do update Triangle::neighbors accordingly.

If user wants to extract boundary vertices, he takes first 3 from first triangle: vertices[0,1,2] , then he takes vertices[0] from every next triangle in the polygon.

My question is, do I need to update some other data inside Triangulation containers, after triangle rotations, so everything else won't get broken?

[artem-ogre]

As long as triangle's winding does not change, you should be fine (I think). But you might want to extract data from Triangulationand discard it (Triangulation) before post-processing it anyway to be 100% sure.

[msokalski]

I keep the original direction (ccw). Thank you very much.

[msokalski]

I've finished 'polygonization' routine. It makes possible to compare results even if there are multiple valid solutions to the same input.
There's only one tiny case where it's not able to 'unify' results. It occurs if we constrain triangulation with crossing edges without adding intersection points. Probably there's nothing we could do to cover it.

Suggestions are welcome!

typedef std::vector<TriInd> Poly;
template <typename T, typename L = LocatorKDTree<T> >
std::vector<Poly> Polygonize(Triangulation<T, L>& cdt)
{
    typedef TriInd PolyInd;
    // for every tri, here we store index of the poly the tri belongs to
    auto map = std::vector<PolyInd>(cdt.triangles.size());
    // vector of Polys, each as a vector of tri indices belonging to the Poly
    auto polys = std::vector<Poly>();

    for (TriInd iT = 0; iT < cdt.triangles.size(); iT++)
    {
        const auto& tri = cdt.triangles[iT];
        auto merged = false;

        // compare i'th tri with its adjacent tris
        for (const auto adj : tri.neighbors)
        {
            // but only if there's adjacent tri and it is processed already
            if (adj == noNeighbor || adj > iT)
                continue;
            // locate reflex vert in adj triangle
            const auto& vr =
                cdt.vertices[opposedVertex(cdt.triangles[adj], iT)];
            const auto& v1 = cdt.vertices[tri.vertices[0]];
            const auto& v2 = cdt.vertices[tri.vertices[1]];
            const auto& v3 = cdt.vertices[tri.vertices[2]];
            using predicates::adaptive::incircle;
            if (!incircle(vr.x, vr.y, v1.x, v1.y, v2.x, v2.y, v3.x, v3.y))
            {
                if (!merged)
                {
                    // append tri to already existing poly
                    merged = true;
                    const auto append_to = map[adj];
                    map[iT] = append_to;
                    polys[append_to].push_back(iT);
                }
                else
                {
                    const auto merge_to = map[iT];
                    const auto merge_from = map[adj];

                    if (merge_to == merge_from)
                        continue;

                    // funny case, tri is a bridge between 2 polys merge'em all
                    // together
                    for (const auto i : polys[merge_from])
                    {
                        map[i] = merge_to;            // remap
                        polys[merge_to].push_back(i); // merge
                    }
                    if (merge_from != polys.size() - 1)
                    {
                        // replace merge_from poly with last poly in polys
                        polys[merge_from] = polys.back();
                        for (const auto i : polys[merge_from])
                        {
                            map[i] = merge_from; // remap
                        }
                    }
                    polys.pop_back();
                }
            }
        }

        if (!merged)
        {
            // at the moment, just alone tri
            // make a new poly for it
            map[iT] = polys.size();
            polys.push_back({ iT });
        }
    }

    // post-proc
    // first triangle in poly should define first 3 vertices at Triangle::vertices[0,1,2]
    // every next triangle should define additional 1 vertex at Triangle::vertices[0]

    // triangle iterator around vertex, 
    // noNeighbor is not checked here!
    struct Stepper
    {
        const Triangulation<T, L>& cdt;
        TriInd current;
        int around;

        Stepper(const Triangulation<T, L>& cdt, TriInd t, int a) : cdt(cdt), current(t), around(a) {}
            
        void StepOver(int a) 
        { 
            around = a; 
        }

        TriInd Clockwise()
        {
            const auto& prev = cdt.triangles[current];

            current = prev.neighbors[around];
            const auto& next = cdt.triangles[current];

            VertInd v = prev.vertices[around];
            if (next.vertices[0] == v)
                around = 0;
            else
            if (next.vertices[1] == v)
                around = 1;
            else
                around = 2;

            return current;
        }
    };

    // bit-mask of a real poly index (without marker)
    const PolyInd mask = (~(PolyInd)0) >> 1;

    for (PolyInd p = 0; p < polys.size(); p++)
    {
        // single triangle polys are ok already
        if (polys[p].size() == 1)
            continue;

        // q = unmarked poly index p
        const PolyInd q = p | ~mask;

        // find good starting triangle,
        // one with exeactly 1 inner edge
        TriInd first = noNeighbor;
        for (const auto t : polys[p])
        {
            if (first == noNeighbor)
            {
                const auto& tri = cdt.triangles[t];
                int inner_edges =
                    (tri.neighbors[0] != noNeighbor && (map[tri.neighbors[0]] & mask) == p) +
                    (tri.neighbors[1] != noNeighbor && (map[tri.neighbors[1]] & mask) == p) +
                    (tri.neighbors[2] != noNeighbor && (map[tri.neighbors[2]] & mask) == p);

                if (inner_edges == 1)
                    first = t;
            }

            // unmark all tris (not inserted to final poly)
            map[t] = q;
        }

        // we can clear current poly now, 
        // as we depend only on map and adjacency
        polys[p].clear();

        TriInd f = first; // current face
        bool step_on = false; // true if current vertex is inserted
        int insert = 2; // 2 for first triangle, all other 0

        Stepper it(cdt, f, 
            cdt.triangles[f].neighbors[0] != noNeighbor && map[cdt.triangles[f].neighbors[0]] == q ? 0 :
            cdt.triangles[f].neighbors[1] != noNeighbor && map[cdt.triangles[f].neighbors[1]] == q ? 1 : 2);

        while (1)
        {
            if (!step_on && map[f] == q)
            {
                step_on = true;
                map[f] = p; // mark as inserted

                if (it.around != insert)
                {
                    // need to rotate
                    auto& tri = cdt.triangles[f];
                    static const int rot[3][3] = { {0,1,2},{2,0,1},{1,2,0} };
                    const int r = rot[it.around][insert];
                    const auto v = tri.vertices[r];
                    const auto n = tri.neighbors[r];
                    switch (r)
                    {
                    case 1:
                        tri.vertices[1] = tri.vertices[0];
                        tri.vertices[0] = tri.vertices[2];
                        tri.vertices[2] = v;
                        tri.neighbors[1] = tri.neighbors[0];
                        tri.neighbors[0] = tri.neighbors[2];
                        tri.neighbors[2] = n;
                        break;
                    case 2:
                        tri.vertices[2] = tri.vertices[0];
                        tri.vertices[0] = tri.vertices[1];
                        tri.vertices[1] = v;
                        tri.neighbors[2] = tri.neighbors[0];
                        tri.neighbors[0] = tri.neighbors[1];
                        tri.neighbors[1] = n;
                        break;
                    default:
                        break;
                    }
                    // adjust iterator after rotation
                    it.StepOver(insert);
                }

                polys[p].push_back(f);

                // everything but first should use 0
                insert = 0;
            }

            TriInd probe = cdt.triangles[f].neighbors[it.around];
            if (probe == noNeighbor || (map[probe] & mask) != p)
            {
                // we're on the last tri inside poly (marked or unmarked)
                assert(step_on); 
                // step on the other leg (ccw in current tri):
                static const int other_leg[3] = { 1,2,0 };
                it.StepOver(other_leg[it.around]);
                step_on = false;
                continue;
            }

            f = it.Clockwise();
            if (f == first)
                break;
        }
    }

    return polys;
}

[msokalski]

Yes I wish to compare results from this gem (CDT) with my own approach called delabella, it's on github (working currently on hybrid branch).

I think it's not a big deal to get particular lib to be consistent, I think CDT even with randomized point insertions produces consistent results.

Now I'm trying to figure out exactly same thing as you mentioned, if it would be possible to define some extra rule to constrained triangulation that when implemented would produce same results across different libs even with input edges crossing.

[artem-ogre]

I see :) Great to meet a fellow library maintainer. This is very interesting! Please let me know if you will have any good findings or ideas (performance, correctness, etc).

I think CDT even with randomized point insertions produces consistent results.
That's because random seed is fixed internally (easier to test and debug this way). But if it wouldn't the results would change from run to run.

Regarding crossing edges: this is explicitly mentioned as not a valid input without using CDT::IntersectingConstraintEdges::Resolve. The output is not guaranteed to be a valid constrained Delaunay triangulation. So I don't think this is anything worth comparing against.

[msokalski]

How about adding compile time option or a template flag, enabling strict triangulation of problematic 'polygons'?
In example, such strict mode would produce triangle fans staring at left-most vertex, if there are 2 we could take lower one.

This may be quite painful for the performance and detecting areas requiring such 'fixes' in the triangulation is probably not trivial.
Also in case of constraining, it should be done after every single edge insertion to prevent previously illustrated case.

EDIT: alternatively we could just pick a vertex with the smallest index if this is consistent with removed duplicates etc.

[artem-ogre]

Here is some interesting related work: https://arxiv.org/pdf/1803.11436.pdf

Given that we detected all polygons formed by concyclic points, we could re-triangulate them with this algorithm that produces a unique solution. This should be done after inserting points and before inserting constraints.

[msokalski]

Thanks, it sounds interesting as aparat from the uniqueness it also improves quality a bit (for use in FE solvers etc.). I'll dig into it on weekend,

[msokalski]

Hi @artem-ogre
So far I didn't find any single (valid) input resulting in 2 different solutions from our 2 libs.
Recently, because my own interval-exact arithmetic thing compared to predicates you use was rather slow I've switched over to it as well.
Currently I'm adding interior / exterior face classifier, so we could also compare on eraseOuterTrianglesAndHoles() results.

[artem-ogre]

Neat! I am away from the keyboard for two weeks. I would be interested in taking a closer look after.

[msokalski]

Would it be possible (reasonably easy) to allow calling eraseOuterTrianglesAndHoles() after eraseSuperTriangle() has been called?
Currently the only option for me is to run CDT 2 times and finalize it with one of erase calls each time.

[artem-ogre]

Would it be possible (reasonably easy) to allow calling eraseOuterTrianglesAndHoles() after eraseSuperTriangle() has been called?
Currently the only option for me is to run CDT 2 times and finalize it with one of erase calls each time.

Triangulation class has value semantics. It can be copied at any point. So you could make a copy after vertices and edges are inserted and then call different erase-methods on two different copies.

Triangulation cdt;
// insert vertices and edges
Triangulation copy = cdt;
cdt.eraseSuperTriangle();
copy.eraseOuterTrianglesAndHoles();

[msokalski]

Ah right, that is too obvious to notice :) Thank you!