Fix issue #25

src/DiscreteVoronoi.jl

@@ -11,7 +11,7 @@ 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

@@ -69,17 +69,28 @@ end
     elseif length(sites) == 1 # Same if there is a single site
         view(grid, TL[1]:BR[1], TL[2]:BR[2]) .= Ref(first(sites))
     else
-        # Otherwise we check if all corners have the same closest site
+        # 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 +104,43 @@ 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))
     else
-        # Otherwise we check if all corners have the same closest site
+        # 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)
+        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 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
+end

src/elimination_methods.jl

@@ -16,23 +16,23 @@ end
 
 function centre_anchor_condition(site, sites, TL, BR; distance=euclidean)
     centre = @. TL + BR / 2
-    anchor = find_closest_site(centre, sites)
+    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
-    anchor = find_closest_site(centre, sites)
+    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 early_stop_sort!(sites, predicate)
+    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 early_stop_sort!(sites, predicate)
-end
+    return unstable_partition!(predicate, sites)[1]
+end

src/helper_functions.jl

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