fast_datashader.jl

"""
    fill_tile!(tile::AbstractMatrix, int_points::AbstractVector{Tuple{Int,Int}}; kw...)

Fill a tile from a region of points converted to `Integer` tuples and sorted.

# Arguments
- `tile`: A square matrix of `Int32`.
- `int_points`: A vector of integer tuple points.

# Keywords
- `x`: x coord of the tile as an Integer.
- `y`: y coord of the tile as an Integer.
- `z`: zoom level
- `zmax`: maximum zoom level
- `tile_size`: size of each tile, default `256`
"""
function fill_tile!(tile::AbstractMatrix, int_points::Vector{Tuple{T,T}}; 
    x=1, y=1, z=0, zmax, tile_size=256
) where T
    zshift = zmax - z
    xmin = T(((x - 1) * tile_size) << zshift + 1)
    ymin = T(((y - 1) * tile_size) << zshift + 1)
    xmax = T((x * tile_size) << zshift)
    ymax = T((y * tile_size) << zshift)
    prev_p_agg = (T(-1), T(-1))
    agg = zero(T)
    i = searchsortedfirst(int_points, (xmin, ymin))
    while i <= lastindex(int_points)
        p = int_points[i]
        p[1] > xmax && break # End of the sorted region
        if p[2] > ymax  
            i = searchsortedfirst(int_points, (p[1] + 1, ymin))
            continue 
        end
        if p[2] < ymin 
            i = searchsortedfirst(int_points, (p[1], ymin))
            continue 
        end
        p_agg = T((p[1] - xmin) >> zshift + oneunit(T)), T((p[2] - ymin) >> zshift + oneunit(T))
        if p_agg === prev_p_agg
            agg += oneunit(T)
        else
            if agg > zero(T)
                tile[prev_p_agg...] += agg
            end
            prev_p_agg = p_agg
            agg = oneunit(T)
        end
        i += 1
    end
    if agg > zero(T)
        tile[prev_p_agg...] += agg
    end
    return tile
end
using ThreadsX

# Generate 100 million points
ps = [(rand(Float32), rand(Float32)) for _ in 1:1e8]

# Scale up, convert to Int32 and Sort 
int_points = @time let
    int_points = Vector{Tuple{Int32,Int32}}(undef, length(ps))
    ThreadsX.map!(int_points, ps) do p
        trunc(Int32, (p[1] * 2^12)) + Int32(1), 
        trunc(Int32, (p[2] * 2^12)) + Int32(1)
    end
    ThreadsX.sort!(int_points) # 4x faster than base sort on 8 cores
end
  # 3.199940 seconds (29.63 k allocations: 2.237 GiB, 0.20% gc time, 4.88% compilation time)

ps = nothing; GC.gc() # Free memory from the original points, if you want to try a billion...

tile = zeros(Int32, 256, 256)
zmax = convert(Int, log2(2^12÷256))
@time fill_tile!(tile, int_points; x=1, y=1, z=0, zmax);
  # 0.312981 seconds (2 allocations: 96 bytes)

# Basic tests for z=0 and z=1
@assert sum(fill_tile!(zeros(Int32, 256, 256), int_points; x=1, y=1, z=0, zmax)) == length(int_points)
let agg = 0
    for x in 1:2, y in 1:2
        agg += sum(fill_tile!(zeros(Int32, 256, 256), int_points; x, y, z=1, zmax))
    end
    @assert agg == length(int_points)
end

using GLMakie
Makie.heatmap(tile)

Ah, I meant something like Extent(X = (0, 1), Y = (0, 1)).

Cool. Just wanted to understand how we correlate zoom levels and tile positions to "real space". I think I get it now but will try to prove it to be sure.