Added points which are on the boundary are added to the constrained segments

[Kevin-Mattheus-Moerman]

@DanielVandH I am working on a minimal example (only have a clunky one currently), but I've noticed that when one has added interior points that happen to be exactly on a boundary edges, they are then included/dded to the constrained segments provided. This seems to me to be undesirable behaviour. It would be better to exclude both "in" and "on" points from the triangulation.

[DanielVandH]

I'm not sure why this is not ideal. The point being added is indeed on a segment and so the segment should get split into two.

This is a very important feature in e.g. mesh refinement to avoid long edges. Similarly in e.g. cell simulations it's not uncommon to have vertices collinear with each other initially, where again this feature is important.

Ultimately if you would like to avoid this feature, it's the users responsibility. I do not intend to change triangulate to ignore user input in this way.

On Wed, 24 Apr 2024, 2:59 pm Kevin Mattheus Moerman, < @.***> wrote:

@DanielVandH https://github.com/DanielVandH I am working on a minimal example (only have a clunky one currently), but I've noticed that when one has added interior points that happen to be exactly on a boundary edges, they are then included/dded to the constrained segments provided. This seems to me to be undesirable behaviour. It would be better to exclude both "in" and "on" points from the triangulation.

— Reply to this email directly, view it on GitHub https://github.com/JuliaGeometry/DelaunayTriangulation.jl/issues/105, or unsubscribe https://github.com/notifications/unsubscribe-auth/AWZPH4HBT4BN4GZQGBXQHVTY663D3AVCNFSM6AAAAABGW7NZQ6VHI2DSMVQWIX3LMV43ASLTON2WKOZSGI3DCMZXG4ZDGMQ . You are receiving this because you were mentioned.Message ID: @.***>

[Kevin-Mattheus-Moerman]

@DanielVandH thanks. I understand what you are saying. However, I think when one uses constraints, I interpret them as "edges" e.g. I want an edge from point n to point m. However if a new point p is added on the edge n<->m then instead we get the edge n<->p and p<->m. I think you what you've implemented is boundary_nodes which form the boundary, rather than boundary edges. So since n and m will be part of the output boundary as per the constrain, you'd say that is fine. Is there a way to perhaps add an option though to say that these are the sole boundary points, i.e. to see them as contrained edges? Perhaps it is a simple optional parameter, e.g. constraint_edges = true? That way it would not alter the behaviour for the other applications you mentioned. I've implemented a fix now (albeit inefficient), so this is not too urgent. But thanks anyway for this great package and your time looking into this.

[Kevin-Mattheus-Moerman]

Below are some test cases for what I built with this in case you are interested:

see also regiontrimesh here.

[DanielVandH]

Hm. I see your point, yes I agree that this could be nice. I think I have an efficient way to implement this. Let me try and get an example working soon and I'll ask for your input.

[DanielVandH]

Haven't been able to get much progress with this since there's a lot of technical details here to do this efficiently (and avoid the annoying amount of edge cases - this package could be so much smaller if not for edge cases :) ), but I did get some ideas down that maybe you are interested in going forward with. I don't really have much time to delve into it much more than this though, sorry.

The idea would be that you first triangulate, and then afterwards you just call enforce_constraints!(tri, segments) (or a keyword argument, as you suggest, would be used to just call enforce_constraints! at the end of triangulate I suppose).

using DelaunayTriangulation: integer_type,
    contains_unoriented_edge,
    is_boundary_node, 
    get_boundary_chain

"""
    fuse_edge!(tri::Triangulation, i, j, r)

Given two vertices `(i, j)` separated only by a single vertex `r`, deletes
that vertex and joins the two vertices with an edge. 
"""
function fuse_edge!(tri::Triangulation, i, j, r) # this is the inverse of split_edge!
    # Problem with this function is that it assumes that the vertex adjoining
    # the edges (i, r) and (r, j) are the same, and similarly for the edges 
    # (j, r) and (r, i). This is not always the case. We need to modify 
    # delete_point! so that it works for boundary points (and for 
    # segment vertices) so that we can just do 
    #   delete_point!(tri, r)
    #   add_segment!(tri, i, j)
    w₁ = get_adjacent(tri, i, r) # assumption: get_adjacent(tri, i, r) == get_adjacent(tri, r, j)
    w₂ = get_adjacent(tri, j, r) # assumption: get_adjacent(tri, j, r) == get_adjacent(tri, r, i)
    delete_triangle!(tri, i, r, w₁)
    delete_triangle!(tri, r, j, w₁)
    delete_triangle!(tri, j, r, w₂)
    delete_triangle!(tri, r, i, w₂)
    add_triangle!(tri, i, j, w₁)
    add_triangle!(tri, j, i, w₂) # need to make sure ghost triangles are handled correctly here 
    return tri
end

"""
    get_segment_chain(tri::Triangulation, i, j) -> Vertices

Given two vertices `i` and `j`, returns the chain of vertices that connects them.
"""
function get_segment_chain(tri::Triangulation, i, j)
    if is_boundary_node(tri, i)[1] && is_boundary_node(tri, j)[1] # that this implies (i, j) should be a boundary edge might not always be true, 
                                                                  # i.e. there might be a segment chain that includes corner points and interior points. 
                                                                  # For the application in this issue, this is fine I believe.
        g = is_boundary_node(tri, i)[2] # the ghost vertex associated with the boundary edge
        chain = get_boundary_chain(tri, i, j, g) # assumes that (i, j) follows the orientation of the boundary - see ?get_boundary_chain
    else 
        chain = get_interior_segment_chain(tri, i, j)
    end     
    return chain
end
function get_interior_segment_chain(tri::Triangulation, i, j)
    I = integer_type(tri)
    chain = [I(i)]
    w = last(chain)
    while w ≠ j # assumes that the chain ends in j. Not always true - e.g. 
                # i could adjoin many segments and the below code will fail
        for k in get_neighbours(tri, w)
            if k ≠ i 
                if contains_segment(tri, w, k) # For the example mentioned above that this could be problematic
                                               # if w adjoins many segments (in which case we may never find j), 
                                               # one way to resolve this for your application would be to check for collinearity, e.g. force 
                                               #     is_collinear(cert),
                                               # where 
                                               #     cert = point_position_relative_to_line(tri, w, k, j)
                                               # If it's possible that the collinearity isn't exact, though, this can 
                                               # still be problematic. You could just try and find all chains until all 
                                               # possible chains end (being careful to not loop the chains), and then 
                                               # take the one that ends in j. Maybe that's what you do in regiontrimesh? I didn't look too much into the details for yours.
                    push!(chain, k)
                    break 
                end
            end # If there are no other segments, note that this will just cause the loop to go on forever.
        end
        w = last(chain)
    end
    return chain
end

"""
    enforce_constraint!(tri::Triangulation, i, j)

Given a segment `(i, j)` represented in `tri` by a sequence of 
collinear segments, removes the vertices defined by the collinear segments 
so that only the segment `(i, j)` remains.
"""
function enforce_constraint!(tri::Triangulation, i, j)
    chain = get_segment_chain(tri, i, j)
    # Note that chain[1] = i and chain[end] = j
    n = length(chain)
    for k in 2:(n-1) # this loop does nothing if contains_unoriented_edge(tri, i, j) is already true, since n = 2
        fuse_edge!(tri, i, chain[k+1], chain[k])
        # Here, note that we originally start with 
        #   i -- chain[2] -- chain[3] -- ... -- chain[end-1] -- j
        # and we want to end with
        #   i -- j
        # So, we start with fusing the chain 
        #   i -- chain[2] -- chain[3]
        # which leaves 
        #   i -- chain[3]
        # or 
        #   i -- chain[3] -- ... -- chain[end-1] -- j
    end
    # Will need to also take some care here, 
    # making sure that all the other segments are removed from 
    #   get_all_segments(tri)
    #   get_interior_segments(tri) (if appropriate)
    # and that (i, j) is added to the appropriate lists.
    return tri 
end

"""
    enforce_constraints!(tri::Triangulation, segments)

Given a list of segments `segments`, enforces the constraints. See [`enforce_constraint`](@ref).
"""
function enforce_constraints!(tri::Triangulation, segments)
    for e in each_edge(segments)
        i, j = edge_vertices(e)
        enforce_constraint!(tri, i, j)
    end
    return tri
end

It would be nice if delete_point! worked for segment vertices for this, but it doesn't at the moment due to there being too many edge cases that I kept encountering I had to give up on it.

[DanielVandH]

Did any of the above work out for you @Kevin-Mattheus-Moerman?

[DanielVandH]

Going to close this due to inactivity and I'm not sure about this feature. The functions I give above maybe work but they won't in general - maybe for your usecase you can take those ideas though.

Feel free to reopen if needed.

[DanielVandH]

I'll re-open this now since maybe I'll try and work on it. This could be added with something like a strict keyword (default false). There are two approaches to this:

1. Reject splitting collinear segments while triangulating.
2. Alternatively, add a post-processing function that fuses edges (like I tried above). I think the first might be the easiest to implement.

The main complication with trying to reject points added onto the boundary is that the triangulations are computed in three stages:

• All the points are added and triangulated (giving the unconstrained Delaunay triangulation).
• Then, all the constrained segments are added.
• Finally, boundaries are enforced so that all exterior points are removed.

Only stages 1-2 are relevant here. The issue is that, since all the points are added in stage 1, it's not as simple as skipping process_collinear_segments in the add_segment! function, since all those points have already been added. If we treat the segment as just being (i, j) without deleting the intersecting points, the triangulation will be invalid. So, process_collinear_segments could take the strict keyword so that any intersected points are deleted. Something like:

function process_collinear_segments!(tri::Triangulation, e, collinear_segments; strict=false, rng::AbstractRNG=Random.default_rng())
    isempty(collinear_segments) && return false
    connect_segments!(collinear_segments)
    extend_segments!(collinear_segments, e)
    if strict
        return _process_collinear_segments_fuse!(tri, collinear_segments; rng)
    else
        return _process_collinear_segments_split!(tri, e, collinear_segments; rng)
    end
end
function _process_collinear_segments_fuse!(tri::Triangulation, collinear_segments; rng::AbstractRNG=Random.default_rng())
    # collinear_segments now looks like
    #   (e[1], i1) - (i1, i2) - (i2, i3) - ⋯ - (in-1, in) - (in, e[2])
    for i in firstindex(collinear_segments):(lastindex(collinear_segments)-1)
        v = terminal(collinear_segments[i])
        delete_point!(tri, v; rng)
    end
    return true
end
function _process_collinear_segments_split!(tri::Triangulation, e, collinear_segments; rng::AbstractRNG=Random.default_rng())
    all_segments = get_all_segments(tri) # the difference between all_segments and segments is important since we need to be careful about what segments are already in the triangulation 
    delete_edge!(all_segments, e)
    split_segment!(tri, e, collinear_segments)
    if contains_boundary_edge(tri, e)
        split_boundary_edge_at_collinear_segments!(tri, collinear_segments)
    end
    for η in each_edge(collinear_segments)
        add_segment!(tri, η; rng)
    end
    return true
end

The only remaining thing that needs to be done to get this working (if this is the right idea, maybe I've overlooked something) is to make delete_point! actually work for points on the boundary. This is not a small fix however.. it did used to work but the edge cases are annoying.

I think this idea is probably the correct one since it involves no post-processing or retriangulating..

[DanielVandH]

I probably can't get around to looking at fixing delete_point! too soon but I'm happy to help go through it if you wanted to have a try @Kevin-Mattheus-Moerman. Otherwise I'll let you know once this might be working.

[DanielVandH]

Another question I'd be interested in is: Why are points being added to your examples that end up falling on the boundary? You will end up with poorer meshes if you try and delete the points. I think it could be worth considering just going forward with these refined edges - all the information you need is in Triangulation already.