JuliaGeometry / DelaunayTriangulation.jl

DelaunayTriangulation.jl: A Julia package for Delaunay triangulations and Voronoi tessellations in the plane
https://juliageometry.github.io/DelaunayTriangulation.jl/
MIT License
70 stars 6 forks source link

Simplify some of the tests, and improve validate_triangulation #51

Closed DanielVandH closed 1 year ago

DanielVandH commented 1 year ago

The tests currently take about an hour to complete. It should be at most 20 minutes I reckon - a lot of time is spent on the same triangulations, using randomised insertion orders to check for issues (it is not impossible that something works 99 times and fails once, I've had it previously). Perhaps we can eliminate this in some places, or somehow use multithreading in some places.

A big issue is the function validate_triangulation in tests/helper_functions.jl, where we check that a constrained triangulation is valid, in particular test_delaunay_criterion:

function test_delaunay_criterion(tri; check_ghost_triangle_delaunay=true)
    for T in each_triangle(tri)
        if DT.is_ghost_triangle(T) && !check_ghost_triangle_delaunay
            continue
        end
        for r in each_solid_vertex(tri)
            cert = DT.point_position_relative_to_circumcircle(tri, T, r)
            if DT.is_inside(cert)
                ace = get_all_constrained_edges(tri)
                if DT.is_empty(ace)
                    flag = !DT.is_inside(cert)
                    if !flag
                        println("Delaunay criterion test failed for the triangle-vertex pair ($T, $r).")
                        return false
                    end
                else # This is extremely slow. Should probably get around to cleaning this up sometime.
                    i, j, k = DT.indices(T)
                    ## We need to see if a subsegment separates T from r. Let us walk from each vertex of T to r.
                    i_history = DT.PointLocationHistory{DT.triangle_type(tri),DT.edge_type(tri),DT.integer_type(tri)}()
                    j_history = DT.PointLocationHistory{DT.triangle_type(tri),DT.edge_type(tri),DT.integer_type(tri)}()
                    k_history = DT.PointLocationHistory{DT.triangle_type(tri),DT.edge_type(tri),DT.integer_type(tri)}()
                    DT.is_boundary_index(i) || jump_and_march(tri, get_point(tri, i); point_indices=nothing, m=nothing, try_points=nothing, k=r, store_history=Val(true), history=i_history)
                    DT.is_boundary_index(j) || jump_and_march(tri, get_point(tri, j); point_indices=nothing, m=nothing, try_points=nothing, k=r, store_history=Val(true), history=j_history)
                    DT.is_boundary_index(k) || jump_and_march(tri, get_point(tri, k); point_indices=nothing, m=nothing, try_points=nothing, k=r, store_history=Val(true), history=k_history)
                    i_walk = DT.is_boundary_index(i) ? Set{DT.triangle_type(tri)}() : i_history.triangles
                    j_walk = DT.is_boundary_index(j) ? Set{DT.triangle_type(tri)}() : j_history.triangles
                    k_walk = DT.is_boundary_index(k) ? Set{DT.triangle_type(tri)}() : k_history.triangles
                    all_edges = Set{NTuple{2,DT.integer_type(tri)}}()
                    E = DT.edge_type(tri)
                    for T in Iterators.flatten((each_triangle(i_walk), each_triangle(j_walk), each_triangle(k_walk)))
                        for e in DT.triangle_edges(T)
                            u, v = DT.edge_indices(e)
                            ee = DT.construct_edge(E, min(u, v), max(u, v))
                            push!(all_edges, (min(u, v), max(u, v)))
                        end
                    end
                    all_constrained_edges = Set{NTuple{2,DT.integer_type(tri)}}()
                    for e in each_edge(ace)
                        u, v = DT.edge_indices(e)
                        ee = DT.construct_edge(E, min(u, v), max(u, v))
                        push!(all_constrained_edges, (min(u, v), max(u, v)))
                    end
                    intersect!(all_edges, all_constrained_edges)
                    flags = [DT.line_segment_intersection_type(tri, initial(e), terminal(e), i, r) for e in each_edge(all_edges) for i in filter(!DT.is_boundary_index, (i, j, k))]
                    flag = !all(DT.is_none, flags) || isempty(flags)
                    if !flag
                        println("Delaunay criterion test failed for the triangle-vertex pair ($T, $r).")
                        return false
                    end
                end
            end
        end
    end
    return true
end

I was rushing at the time that I wrote this. I reckon we need only checking the surrounding_polygon for the points, checking if the line from a point to another intersects the polygon at all, and then testing more carefully inside the polygon. triangle_line_segment_intersection would be good for this also - that's a big reason why I initially wrote that predicate...

DanielVandH commented 1 year ago

For reference, here's an example test run:

Test Summary:                   |      Pass      Total      Time
DelaunayTriangulation           | 157101159  157101159  44m19.0s
  Triangulation                 |    300728     300728  23m34.2s
    Gmsh                        |     61005      61005     55.6s
    Rectangular Triangulation   |         3          3      4.0s
    Bowyer Watson               |      2123       2123     32.2s
    Triangulate                 |      2302       2302     37.4s
    Convex Triangulation        |    174940     174940      8.9s
    Constrained Triangulation   |     60355      60355  21m16.0s
  Utilities                     |     51272      51272     33.1s
    Utilities                   |     40038      40038      5.5s
    Geometric Utilities         |     11206      11206      2.8s
    Triangulation Validation    |        28         28     24.8s
  Interfaces                    |      1502       1502      3.5s
    Triangles                   |       151        151      1.5s
    Edges                       |       122        122      1.1s
    Points                      |      1097       1097      0.4s
    Boundary Nodes              |       132        132      0.6s
  Data Structures               |   4815979    4815979     30.5s
    Adjacent                    |       186        186      0.9s
    Adjacent2Vertex             |       190        190      1.1s
    Graph                       |        65         65      0.2s
    ConvexHull                  |     15012      15012      1.6s
    Triangulation               |   4058275    4058275     19.3s
    Representative Points       |       104        104      0.7s
    Statistics                  |    742147     742147      6.6s
  Operations                    |   8333983    8333983   4m46.8s
    add_triangle!               |        79         79      2.4s
    delete_triangle!            |       538        538      1.1s
    add_ghost_triangles!        |        72         72      1.0s
    delete_ghost_triangles!     |     53918      53918      4.6s
    add_point!                  |       503        503     11.0s
    flip_edge!                  |        13         13      0.3s
    split_triangle!             |   1250012    1250012     11.5s
    split_edge!                 |   7002408    7002408     38.0s
    legalise_edge!              |        35         35      0.3s
    delete_point!               |     14587      14587   3m28.5s
    (un)lock_convex_hull!.jl    |       912        912      2.1s
    delete_holes!               |     10906      10906      6.2s
  Predicates                    |    578880     578880     15.3s
    Certificate                 |       424        424      0.1s
    Boundaries and Ghosts       |       984        984      1.3s
    General                     |    572942     572942      6.2s
    Index and Ghost Handling    |      4530       4530      7.8s
  Point Location                | 124868197  124868197   3m49.5s
    Brute Force                 |       160        160      0.4s
    Initial Point Selection     |       185        185      0.8s
    Initial Triangle Selection  |      4774       4774      0.8s
    Interior Edge Intersections |        74         74      1.0s
    Ghost Search                |       248        248      0.6s
    Jump and March              | 124862756  124862756   3m45.8s
  Constrained Triangulation     |    135894     135894     17.6s
    Segment Location            |    135876     135876     16.4s
    Segment Insertion           |        18         18      1.2s
  Refinement                    |  17422327   17422327   7m38.4s
    Encroachment                |    259618     259618     58.5s
    Quality Assessment          |       141        141      1.7s
    Refinement Operations       |     42009      42009      6.7s
    Refinement                  |  17120559   17120559   6m31.4s
  Documentation images          |    182652     182652     56.3s
  Voronoi                       |    409745     409745   1m53.8s
DanielVandH commented 1 year ago

Some improvements to this are given in #56, which makes validate_triangulation use @threads.

DanielVandH commented 1 year ago
DelaunayTriangulation           | 99449013  99449013  21m12.1s
  Triangulation                 |   298907    298907  10m46.3s
    Gmsh                        |    61005     61005     23.8s
    Rectangular Triangulation   |        3         3      4.2s
    Bowyer Watson               |     2123      2123     15.6s
    Triangulate                 |     2302      2302     21.0s
    Convex Triangulation        |   175253    175253      8.0s
    Constrained Triangulation   |    58221     58221   9m33.6s
  Utilities                     |    51272     51272     23.7s
    Utilities                   |    40038     40038      7.2s
    Geometric Utilities         |    11206     11206      3.0s
    Triangulation Validation    |       28        28     13.5s
  Interfaces                    |     1502      1502      5.0s
    Triangles                   |      151       151      2.0s
    Edges                       |      122       122      1.4s
    Points                      |     1097      1097      0.7s
    Boundary Nodes              |      132       132      0.9s
  Data Structures               |  4496540   4496540     34.8s
    Adjacent                    |      148       148      1.2s
    Adjacent2Vertex             |      190       190      1.3s
    Graph                       |       65        65      0.3s
    ConvexHull                  |    15012     15012      2.0s
    Triangulation               |  3735094   3735094     21.4s
    Representative Points       |      104       104      0.8s
    Statistics                  |   745927    745927      7.6s
  Operations                    | 11406427  11406427   2m01.3s
    add_triangle!               |       79        79      2.7s
    delete_triangle!            |      538       538      1.2s
    add_ghost_triangles!        |       72        72      1.0s
    delete_ghost_triangles!     |    53918     53918      4.8s
    add_point!                  |  3072947   3072947     14.5s
    flip_edge!                  |       13        13      0.3s
    split_triangle!             |  1250012   1250012      6.9s
    split_edge!                 |  7002408   7002408     26.1s
    legalise_edge!              |       35        35      0.3s
    delete_point!               |    14587     14587     55.4s
    (un)lock_convex_hull!.jl    |      912       912      1.3s
    delete_holes!               |    10906     10906      6.8s
  Predicates                    |   578881    578881     17.2s
    Certificate                 |      424       424      0.1s
    Boundaries and Ghosts       |      984       984      1.4s
    General                     |   572942    572942      7.5s
    Index and Ghost Handling    |     4531      4531      8.2s
  Point Location                | 62463205  62463205   1m23.5s
    Brute Force                 |      160       160      0.4s
    Initial Point Selection     |      185       185      2.3s
    Initial Triangle Selection  |     4774      4774      0.9s
    Interior Edge Intersections |       74        74      1.1s
    Ghost Search                |      248       248      0.7s
    Jump and March              | 62457764  62457764   1m18.1s
  Constrained Triangulation     |    70894     70894     10.5s
    Segment Location            |    70876     70876      9.6s
    Segment Insertion           |       18        18      0.9s
  Refinement                    | 19488448  19488448   2m18.0s
    Encroachment                |   259618    259618      1.8s
    Quality Assessment          |      141       141      1.7s
    Refinement Operations       |    42009     42009      6.3s
    Refinement                  | 19186680  19186680   2m08.1s
  Documentation images          |   183192    183192     52.5s
  Voronoi                       |   409745    409745   2m19.2s
     Testing DelaunayTriangulation tests passed