Fix issue #25

benchmark/btime.jl

@@ -0,0 +1,61 @@
+using DiscreteVoronoi
+using BenchmarkTools
+using Random
+using StaticArrays
+
+using DiscreteVoronoi: Coord, exact_aux, centre_anchor_aux, exact_aux_es, centre_anchor_aux_es
+
+function random_coordinates(N, M, K)
+    Coord.(shuffle!([Iterators.product(1:N, 1:M)...])[1:min(N * M, K)])
+end
+
+for n in [400]
+    for s in [isqrt(n), n, n * isqrt(n), n * n]
+        @show n, s
+
+        println("jfa_voronoi!")
+        @btime jfa_voronoi!(grid, sites) setup = (
+            Random.seed!(42);
+            grid = zeros(Coord, ($n, $n));
+            sites = random_coordinates(size(grid)..., $s)) evals=1
+
+        if s <= n
+            println("dac_voronoi!")
+            @btime dac_voronoi!(grid, sites) setup=(
+                Random.seed!(42);
+                grid = zeros(Coord, ($n, $n));
+                sites = random_coordinates(size(grid)..., $s);
+                preset_voronoi!(grid, sites)) evals=1
+        end
+
+        println("redac_voronoi!")
+        @show exact_aux
+        @btime redac_voronoi!(grid, sites, auxiliary=exact_aux) setup=(
+            Random.seed!(42);
+            grid = zeros(Coord, ($n, $n));
+            sites = random_coordinates(size(grid)..., $s);
+            preset_voronoi!(grid, sites)) evals=1
+
+        @show centre_anchor_aux
+        @btime redac_voronoi!(grid, sites, auxiliary=centre_anchor_aux) setup=(
+            Random.seed!(42);
+            grid = zeros(Coord, ($n, $n));
+            sites = random_coordinates(size(grid)..., $s);
+            preset_voronoi!(grid, sites)) evals=1
+
+        println("redac_voronoi_es!")
+        @show exact_aux_es
+        @btime redac_voronoi_es!(grid, sites, auxiliary=exact_aux_es) setup=(
+            Random.seed!(42);
+            grid = zeros(Coord, ($n, $n));
+            sites = random_coordinates(size(grid)..., $s);
+            preset_voronoi!(grid, sites)) evals=1
+
+        @show centre_anchor_aux_es
+        @btime redac_voronoi_es!(grid, sites, auxiliary=centre_anchor_aux_es) setup=(
+            Random.seed!(42);
+            grid = zeros(Coord, ($n, $n));
+            sites = random_coordinates(size(grid)..., $s);
+            preset_voronoi!(grid, sites)) evals=1
+    end
+end

src/DiscreteVoronoi.jl

@@ -3,15 +3,14 @@ module DiscreteVoronoi
 #TODO: Implement hybrid algorithms for certain gridsize and number of sites.
 #TODO: Add nice UnicodePlot recipe :)
 
-export find_closest_site, find_closest_site!, get_corners, get_quadrants, label_voronoi_grid, voronoi_equality # Helper functions
-export naive_voronoi!, jfa_voronoi!, dac_voronoi!, redac_voronoi! # Core functionality
+export unstable_partition!, find_closest_site, find_closest_site!, get_corners, get_quadrants, label_voronoi_grid, voronoi_equality # Helper functions
+export naive_voronoi!, jfa_voronoi!, preset_voronoi!, dac_voronoi!, redac_voronoi!, redac_voronoi_es! # Core functionality
 
-using LinearAlgebra: norm
 using Random: shuffle
 using StaticArrays
 using Distances: euclidean
 
-Coord = SVector{2,Int}
+const Coord = SVector{2,Int}
 
 # Completely non-exported files with core definitions
 include("EarlyStopper.jl") # Because every package needs at least one struct...

src/core_algorithms.jl

@@ -50,6 +50,15 @@ function jfa_voronoi!(grid::Matrix{T}, sites::Vector{T}; distance=euclidean) whe
     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}}
 
@@ -67,19 +76,31 @@ end
     if all(BR .== TL)
         grid[TL...] = find_closest_site(TL, sites, distance=distance)
     elseif length(sites) == 1 # Same if there is a single site
-        view(grid, TL[1]:BR[1], TL[2]:BR[2]) .= Ref(first(sites))
+        grid[TL[1]:BR[1], TL[2]:BR[2]] .= Ref(first(sites))
     else
-        # Otherwise we check if all corners have the same closest site
+        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)
-        closest_corners = (find_closest_site(corner, sites, distance=distance) for corner in corners)
-        if allequal(closest_corners)
-            grid[TL[1]:BR[1], TL[2]:BR[2]] .= Ref(first(closest_corners))
-        else
-            # 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)
+        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
@@ -93,32 +114,58 @@ 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, EarlyStopper(sites); distance=distance, auxiliary=auxiliary)
+    _redac_voronoi!(grid, TL, BR, sites; distance=distance, auxiliary=auxiliary)
     return nothing
 end
 
-@inbounds function _redac_voronoi!(grid, TL, BR, sites::ES; distance, auxiliary) where {ES<:EarlyStopper}
+@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
-        view(grid, TL[1]:BR[1], TL[2]:BR[2]) .= Ref(first(sites))
+        grid[TL[1]:BR[1], TL[2]:BR[2]] .= Ref(first(sites))
     else
-        # Otherwise we check if all corners have the same closest site
-        corners = get_corners(TL, BR)
-        closest_corners = (find_closest_site!(grid, corner, sites, distance=distance) for corner in corners)
-        if allequal(closest_corners)
-            view(grid, TL[1]:BR[1], TL[2]:BR[2]) .= Ref(first(closest_corners))
-        else
-            # 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
+        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
+end

src/elimination_methods.jl

@@ -15,14 +15,30 @@ function exact_condition(site, sites, TL, BR; distance=euclidean)
 end
 
 function centre_anchor_condition(site, sites, TL, BR; distance=euclidean)
-    centre = @. TL + BR / 2
-    anchor = find_closest_site(centre, sites)
+    centre = @. (TL + BR) / 2
+    anchor = find_closest_site(centre, sites; distance=distance)
     corners = get_corners(TL, BR)
     return any(distance(site, corner) <= distance(anchor, corner) for corner in corners)
 end
 
 function centre_anchor_aux(grid, sites, TL, BR; distance=euclidean)
-    centre = @. TL + BR / 2
+    centre = @. (TL + BR) / 2
+    anchor = find_closest_site(centre, sites; distance=distance)
+    corners = get_corners(TL, BR)
+    function predicate(site)
+        return all(TL .<= site .<= BR) ||
+               any(distance(corner, site) <= distance(corner, anchor) for corner in corners)
+    end
+    return unstable_partition!(predicate, sites)[1]
+end
+
+function exact_aux(grid, sites, TL, BR; distance=euclidean)
+    predicate(site) = all(TL .<= site .<= BR) || exact_condition(site, sites, TL, BR; distance=distance)
+    return unstable_partition!(predicate, sites)[1]
+end
+
+function centre_anchor_aux_es(grid, sites, TL, BR; distance=euclidean)
+    centre = @. (TL + BR) / 2
     anchor = find_closest_site(centre, sites)
     corners = get_corners(TL, BR)
     function predicate(site)
@@ -32,7 +48,7 @@ function centre_anchor_aux(grid, sites, TL, BR; distance=euclidean)
     return early_stop_sort!(sites, predicate)
 end
 
-function exact_aux(grid, sites, TL, BR; distance=euclidean)
+function exact_aux_es(grid, sites, TL, BR; distance=euclidean)
     predicate(site) = all(TL .<= site .<= BR) || exact_condition(site, sites, TL, BR; distance=distance)
     return early_stop_sort!(sites, predicate)
-end
+end

src/helper_functions.jl

@@ -22,6 +22,27 @@ function swap!(v, i::Int, j::Int)
     return nothing
 end
 
+@inbounds function unstable_partition!(f, A)
+    i, j = 1, length(A)
+    while i <= j
+        while i <= j && f(A[i])
+            i += 1
+        end
+        # @assert i > length(A) || !f(A[i])
+        i >= j && @views return A[1:(i-1)], A[i:length(A)]
+        while i < j && !f(A[j])
+            j -= 1
+        end
+        # @assert i == j || f(A[j])
+        i == j && @views return A[1:(i-1)], A[i:length(A)]
+        A[i], A[j] = A[j], A[i]
+        i += 1
+        j -= 1
+    end
+    # @assert i == j + 1
+    @views return A[1:j], A[i:length(A)]
+end
+
 """
     get_quadrants(TL, BR)
 

test/runtests.jl

@@ -7,7 +7,8 @@ using Distances
 
 include("helper_functions.jl")
 
-using DiscreteVoronoi: get_corners, voronoi_equality, exact_condition, exact_aux, centre_anchor_aux
+using DiscreteVoronoi: get_corners, voronoi_equality, exact_condition
+using DiscreteVoronoi: exact_aux, centre_anchor_aux, exact_aux_es, centre_anchor_aux_es
 using DiscreteVoronoi: EarlyStopper, early_stop_sort!
 
 const Coord = SVector{2,Int}
@@ -120,6 +121,35 @@ end
 
 const REPETITIONS = 100
 
+@testset "unstable_partition! working" begin
+    function test(site, N, M) 
+        global counter
+        counter += 1
+        (site[1] < N ÷ 2 && site[2] < M ÷ 2) || 
+        (site[1] > N ÷ 2 && site[2] > M ÷ 2)
+    end
+    Random.seed!(42)
+    test_passes = true
+    for i in 1:REPETITIONS
+        N, M = rand(1:100, 2)
+        sites = random_coordinates(N, M, rand(1:100))
+
+        global counter = 0
+        trues, falses = unstable_partition!(sites) do site
+            test(site, N, M)
+        end
+        test_passes = test_passes && counter == length(sites)
+        for site in sites
+            if test(site, N, M) 
+                test_passes = test_passes && site in trues 
+            else 
+                test_passes = test_passes && site in falses
+            end
+        end
+    end
+    test_passes ? (@test test_passes) : (@test_broken test_passes)
+end
+
 @testset verbose=true "jfa_voronoi! working for Int sites" begin
     Random.seed!(42)
     for distance in [cityblock, euclidean, chebyshev]
@@ -186,3 +216,28 @@ end
         end
     end
 end
+
+@testset verbose=true "redac_voronoi_es! matching results for Int sites" begin
+    Random.seed!(42)
+    for distance in [cityblock, euclidean, chebyshev]
+        @testset verbose=true "$distance" begin
+            for auxiliary in [exact_aux_es, centre_anchor_aux_es]
+                @testset "$auxiliary" begin
+                    test_passes = true
+                    for i in 1:REPETITIONS
+                        N, M = rand(1:100, 2)
+                        points = random_coordinates(N, M, rand(1:100))
+
+                        grid1 = zeros(Coord, (N, M))
+                        naive_voronoi!(grid1, points, distance=distance)
+                        grid2 = zeros(Coord, (N, M))
+                        redac_voronoi_es!(grid2, points, distance=distance, auxiliary=auxiliary)
+                        test_passes = test_passes && voronoi_equality(grid2, grid1; distance=distance)
+                        test_passes || break
+                    end
+                    test_passes ? (@test test_passes) : (@test_broken test_passes)
+                end
+            end
+        end
+    end
+end