DiscreteVoronoi.jl/src/core_algorithms.jl at 0c0e3d8322535980c139466441dcae80ae81a9d5 · jacobusmmsmit/DiscreteVoronoi.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
# This is the core functionality for the methods of producing discrete Voronoi
# diagrams. The JFA is most commonly used on a GPU where every cell is run in
# parallel, but it is known to produce (small) errors. The most efficient
# algorithm for single-core performance is redac_voronoi!, but it cannot be
# multithreaded in its current state as it heavily mutates its arguments.

"""
    naive_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean) where {T<:SVector{2,Int}}

Construct in-place a Voronoi diagram in the most basic way possible: check every cell and every combination.
"""
function naive_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean) where {T<:Coord}
    for I in CartesianIndices(grid)
        @inbounds grid[I] = find_closest_site(Tuple(I), sites; distance=distance)
    end
    return nothing
end

"""
    jfa_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean) where {T<:SVector{2,Int}}

Construct in-place a Voronoi diagram using the [jump flooding algorithm](https://en.wikipedia.org/wiki/Jump_flooding_algorithm).
The algorithm assumes that a blank cell in the grid has value `SVector(0, 0)` and that sites are inside the grid.
"""
function jfa_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean) where {T<:Coord}
    for site in sites
        # Splatting (grid[site...] = site) causes allocations?
        x, y = site
        grid[x, y] = site
    end
    k = max(size(grid)...)
    while k > 1
        k = k ÷ 2 + k % 2
        @inbounds for I in CartesianIndices(grid)
            x, y = Tuple(I)
            for j in (-k, 0, k), i in (-k, 0, k)
                checkbounds(Bool, grid, x + i, y + j) || continue
                i == j == 0 && continue
                siteq = grid[x+i, y+j]
                siteq !== SVector(0, 0) || continue
                sitep = grid[x, y]
                if sitep == SVector(0, 0)
                    grid[x, y] = siteq
                elseif distance(sitep, (x, y)) > distance(siteq, (x, y))
                    grid[x, y] = siteq
                end
            end
        end
    end
    return nothing
end

@inbounds function preset_voronoi!(grid, sites)
    for site in sites
        if checkbounds(Bool, grid, site...)
            grid[site...] = site
        end
    end
    return nothing
end

"""
    dac_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean) where {T<:SVector{2,Int}}

Construct in-place a Voronoi diagram using the [divide-and-conquer algorithm](https://doi.org/10.1109%2Feit48999.2020.9208270).
"""
function dac_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean) where {T<:Coord}
    TL = (1, 1)
    BR = size(grid)
    _dac_voronoi!(grid, TL, BR, sites, distance=distance)
    return nothing
end

@inbounds function _dac_voronoi!(grid, TL, BR, sites; distance)
    any(TL .> BR) && return nothing
    if all(BR .== TL)
        grid[TL...] = find_closest_site(TL, sites, distance=distance)
    elseif length(sites) == 1 # Same if there is a single site
        grid[TL[1]:BR[1], TL[2]:BR[2]] .= Ref(first(sites))
    else
        all(site -> site != zero(Coord), @view grid[TL[1]:BR[1], TL[2]:BR[2]]) && return nothing
        # Otherwise we check if all corners have the same closest site ...
        corners = get_corners(TL, BR)
        mins = ((findmin(sites) do site
            distance(corner, site)
        end for corner in corners)...,)
        if allequal(site for (dist, site) in mins)
            # ... and if the closest site is unique
            min_site = sites[mins[1][2]]
            dists = ((minimum(site for site in sites if site != min_site) do site
                distance(corner, site)
            end for corner in corners)...,)
            if all(zip(mins, dists)) do ((min_dist, _), dist)
                    dist > min_dist
                end
                grid[TL[1]:BR[1], TL[2]:BR[2]] .= Ref(min_site)
                return nothing
            end
        end
        # And if not we divide the grid into quadrants and "conquer" each one
        for (quadrant_TL, quadrant_BR) in get_quadrants(TL, BR)
            _dac_voronoi!(grid, quadrant_TL, quadrant_BR, sites; distance=distance)
        end
    end
    return nothing
end

"""
    redac_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean, auxiliary=exact_aux) where {T<:SVector{2,Int}}

Performs a divide-and-conquer method similar to `dac_voronoi!` but has an additional site-elimination
step which aims to reduce the work of subsequent steps.
"""
function redac_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean, auxiliary=exact_aux) where {T<:Coord}
    TL = 1, 1
    BR = size(grid)
    _redac_voronoi!(grid, TL, BR, sites; distance=distance, auxiliary=auxiliary)
    return nothing
end

@inbounds function _redac_voronoi!(grid, TL, BR, sites; distance, auxiliary)
    any(TL .> BR) && return nothing
    # Then, if the grid is a single cell then we are done
    if all(BR .== TL)
        find_closest_site!(grid, TL, sites; distance=distance)
    elseif length(sites) == 1 # Same if there is a single site
        grid[TL[1]:BR[1], TL[2]:BR[2]] .= Ref(first(sites))
    else
        all(site -> site != zero(Coord), @view grid[TL[1]:BR[1], TL[2]:BR[2]]) && return nothing
        # And if not we eliminate faraway seeds from subsequent steps
        # `auxiliary` sorts sites by whether the predicate is true and stores how many are true.
        local_sites = auxiliary(grid, sites, TL, BR; distance=distance)
        # then divide the grid into quadrants and "conquer" each one
        for (quadrant_TL, quadrant_BR) in get_quadrants(TL, BR)
            _redac_voronoi!(grid, quadrant_TL, quadrant_BR, local_sites; distance=distance, auxiliary=auxiliary)
        end
    end
    return nothing
end

"""
    redac_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean, auxiliary=exact_aux) where {T<:SVector{2,Int}}
Performs a divide-and-conquer method similar to `dac_voronoi!` but has an additional site-elimination
step which aims to reduce the work of subsequent steps.
"""
function redac_voronoi_es!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean, auxiliary=exact_aux) where {T<:Coord}
    TL = 1, 1
    BR = size(grid)
    _redac_voronoi_es!(grid, TL, BR, EarlyStopper(sites); distance=distance, auxiliary=auxiliary)
    return nothing
end

@inbounds function _redac_voronoi_es!(grid, TL, BR, sites::ES; distance, auxiliary) where {ES<:EarlyStopper}
    any(TL .> BR) && return nothing
    # Then, if the grid is a single cell then we are done
    if all(BR .== TL)
        find_closest_site!(grid, TL, sites; distance=distance)
    elseif length(sites) == 1 # Same if there is a single site
        grid[TL[1]:BR[1], TL[2]:BR[2]] .= Ref(first(sites))
    else
        all(site -> site != zero(Coord), @view grid[TL[1]:BR[1], TL[2]:BR[2]]) && return nothing
        # And if not we eliminate faraway seeds from subsequent steps
        # `auxiliary` sorts sites by whether the predicate is true and stores how many are true.
        local_sites = auxiliary(grid, sites, TL, BR; distance=distance)
        # then divide the grid into quadrants and "conquer" each one
        for (quadrant_TL, quadrant_BR) in get_quadrants(TL, BR)
            _redac_voronoi_es!(grid, quadrant_TL, quadrant_BR, local_sites; distance=distance, auxiliary=auxiliary)
        end
    end
    return nothing
end