diff --git a/src/screen.jl b/src/screen.jl
index fee8493..a401bee 100644
--- a/src/screen.jl
+++ b/src/screen.jl
@@ -410,6 +410,38 @@ function global_gl_screen(resolution::Tuple, visibility::Bool, tries = 1)
     screen
 end
 
+
+
+#################################################################################
+### Point picking
+################################################################################
+
+
+
+function pick_native(screen::Screen, rect::IRect2D)
+    isopen(screen) || return Matrix{SelectionID{Int}}(undef, 0, 0)
+    window_size = widths(screen)
+    fb = screen.framebuffer
+    buff = fb.buffers[:objectid]
+    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1])
+    glReadBuffer(GL_COLOR_ATTACHMENT1)
+    rx, ry = minimum(rect)
+    rw, rh = widths(rect)
+    w, h = window_size
+    sid = zeros(SelectionID{UInt32}, widths(rect)...)
+    if rx > 0 && ry > 0 && rx + rw <= w && ry + rh <= h
+        glReadPixels(rx, ry, rw, rh, buff.format, buff.pixeltype, sid)
+        for i in eachindex(sid)
+            if sid[i][2] > 0x3f800000
+                sid[i] = SelectionID(0, sid[i].index)
+            end
+        end
+        return sid
+    else
+        error("Pick region $rect out of screen bounds ($w, $h).")
+    end
+end
+
 function pick_native(screen::Screen, xy::Vec{2, Float64})
     isopen(screen) || return SelectionID{Int}(0, 0)
     sid = Base.RefValue{SelectionID{UInt32}}()
@@ -424,80 +456,93 @@ function pick_native(screen::Screen, xy::Vec{2, Float64})
         glReadPixels(x, y, 1, 1, buff.format, buff.pixeltype, sid)
         return convert(SelectionID{Int}, sid[])
     end
-    return return SelectionID{Int}(0, 0)
+    return SelectionID{Int}(0, 0)
 end
 
-function pick_native(screen::Screen, xy::Vec{2, Float64}, range::Float64)
-    isopen(screen) || return SelectionID{Int}(0, 0)
-    window_size = widths(screen)
-    w, h = window_size
-    if !((1.0 <= xy[1] <= w) && (1.0 <= xy[2] <= h))
-        return SelectionID{Int}(0, 0)
+function AbstractPlotting.pick(scene::SceneLike, screen::Screen, xy::Vec{2, Float64})
+    sid = pick_native(screen, xy)
+    if haskey(screen.cache2plot, sid.id)
+        plot = screen.cache2plot[sid.id]
+        return (plot, sid.index)
+    else
+        return (nothing, 0)
     end
+end
 
-    fb = screen.framebuffer
-    buff = fb.buffers[:objectid]
-    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1])
-    glReadBuffer(GL_COLOR_ATTACHMENT1)
+function AbstractPlotting.pick(scene::SceneLike, screen::Screen, rect::IRect2D)
+    map(pick_native(screen, rect)) do sid
+        if haskey(screen.cache2plot, sid.id)
+            (screen.cache2plot[sid.id], sid.index)
+        else
+            (nothing, sid.index)
+        end
+    end
+end
+
+
+# Skips one set of allocations
+function AbstractPlotting.pick_closest(scene::SceneLike, screen::Screen, xy, range)
+    isopen(screen) || return (nothing, 0)
+    w, h = widths(screen)
+    ((1.0 <= xy[1] <= w) && (1.0 <= xy[2] <= h)) || return (nothing, 0)
 
     x0, y0 = max.(1, floor.(Int, xy .- range))
     x1, y1 = min.([w, h], floor.(Int, xy .+ range))
     dx = x1 - x0; dy = y1 - y0
-    sid = Matrix{SelectionID{UInt32}}(undef, dx, dy)
-    glReadPixels(x0, y0, dx, dy, buff.format, buff.pixeltype, sid)
+    sid = pick_native(screen, IRect2D(x0, y0, dx, dy))
 
     min_dist = range^2
     id = SelectionID{Int}(0, 0)
     x, y =  xy .+ 1 .- Vec2f0(x0, y0)
     for i in 1:dx, j in 1:dy
         d = (x-i)^2 + (y-j)^2
-        if (d < min_dist) && (sid[i, j][1] > 0x00000000) && (sid[i, j][2] < 0x3f800000)
+        if (d < min_dist) && (sid[i, j][1] > 0x00000000) && 
+            (sid[i, j][2] < 0x3f800000) && haskey(screen.cache2plot, sid[i, j][1])
             min_dist = d
             id = convert(SelectionID{Int}, sid[i, j])
         end
     end
-    return id
-end
 
-function AbstractPlotting.pick(scene::SceneLike, screen::Screen, xy::Vec{2, Float64})
-    sid = pick_native(screen, xy)
-    if haskey(screen.cache2plot, sid.id)
-        plot = screen.cache2plot[sid.id]
-        return (plot, sid.index)
+    if haskey(screen.cache2plot, id[1])
+        return (screen.cache2plot[id[1]], id[2])
     else
         return (nothing, 0)
     end
 end
 
-function AbstractPlotting.pick(scene::SceneLike, screen::Screen, xy::Vec{2, Float64}, range::Float64)
-    sid = pick_native(screen, xy, range)
-    if haskey(screen.cache2plot, sid.id)
-        plot = screen.cache2plot[sid.id]
-        return (plot, sid.index)
-    else
-        return (nothing, 0)
+# Skips some allocations
+function AbstractPlotting.pick_sorted(scene::SceneLike, screen::Screen, xy, range)
+    isopen(screen) || return (nothing, 0)
+    w, h = widths(screen)
+    if !((1.0 <= xy[1] <= w) && (1.0 <= xy[2] <= h))
+        return Tuple{AbstractPlot, Int}[]
     end
-end
+    x0, y0 = max.(1, floor.(Int, xy .- range))
+    x1, y1 = min.([w, h], floor.(Int, xy .+ range))
+    dx = x1 - x0; dy = y1 - y0
+    
+    picks = pick_native(screen, IRect2D(x0, y0, dx, dy))
 
-function AbstractPlotting.pick(screen::Screen, rect::IRect2D)
-    window_size = widths(screen)
-    fb = screen.framebuffer
-    buff = fb.buffers[:objectid]
-    glBindFramebuffer(GL_FRAMEBUFFER, fb.id[1])
-    glReadBuffer(GL_COLOR_ATTACHMENT1)
-    x, y = minimum(rect)
-    rw, rh = widths(rect)
-    w, h = window_size
-    sid = zeros(SelectionID{UInt32}, widths(rect)...)
-    if x > 0 && y > 0 && x <= w && y <= h
-        glReadPixels(x, y, rw, rh, buff.format, buff.pixeltype, sid)
-        sid = filter(x -> x.id < 0x3f800000,sid)
-        return map(unique(vec(SelectionID{Int}.(sid)))) do sid
-            (screen.cache2plot[sid.id], sid.index)
+    selected = filter(x -> x[1] > 0 && haskey(screen.cache2plot, x[1]), unique(vec(picks)))
+    distances = [range^2 for _ in selected]
+    x, y =  xy .+ 1 .- Vec2f0(x0, y0)
+    for i in 1:dx, j in 1:dy
+        if picks[i, j][1] > 0
+            d = (x-i)^2 + (y-j)^2
+            i = findfirst(isequal(picks[i, j]), selected)
+            if i === nothing
+                @warn "This shouldn't happen..."
+            elseif distances[i] > d
+                distances[i] = d
+            end
         end
     end
-    return Tuple{AbstractPlot, Int}[]
+
+    idxs = sortperm(distances)
+    permute!(selected, idxs)
+    return map(id -> (screen.cache2plot[id[1]], id[2]), selected)
 end
 
+
 pollevents(::GLScreen) = nothing
 pollevents(::Screen) = GLFW.PollEvents()
diff --git a/test/unit_tests.jl b/test/unit_tests.jl
index 6e04515..1e0e800 100644
--- a/test/unit_tests.jl
+++ b/test/unit_tests.jl
@@ -42,13 +42,12 @@ end
         origin_px = project_sp(ax.scene, Point(origin(rect)))
         tip_px = project_sp(ax.scene, Point(origin(rect) .+ widths(rect)))
         rect_px = IRect2D(round.(origin_px), round.(tip_px .- origin_px))
-        #! there is no pick(::Scene,::IRect2D)
-        plot_idx = pick(screen, rect_px)
+        picks = unique(pick(ax.scene, rect_px))
 
         # objects returned in plot_idx should be either grid lines (i.e. LineSegments) or Scatter points
-        @test all(pi-> pi[1] isa Union{LineSegments,Scatter, AbstractPlotting.Mesh}, plot_idx)
+        @test all(pi-> pi[1] isa Union{LineSegments,Scatter, AbstractPlotting.Mesh}, picks)
         # scatter points should have indices equal to those in 99991:99998
-        scatter_plot_idx = filter(pi -> pi[1] isa Scatter, plot_idx)
+        scatter_plot_idx = filter(pi -> pi[1] isa Scatter, picks)
         @test Set(last.(scatter_plot_idx)) == Set(99991:99998)
     end
 end