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)

Just to clarify...this currently operates on ((0, 1), (0, 1)) for tile (1, 1)?
That means that we need to actually figure out how to link "zoom" levels with data limits - or is that something which Tyler might do?

FWIW, negative z values seem to work fine. So the capability is there - we just need to figure out the mechanism.

Alternatively - we could rescale the data at the point we cache it.

Just to clarify...this currently operates on ((0, 1), (0, 1)) for tile (1, 1)?

Nor sure I totally understand this, what is ((0, 1), (0, 1)) ?

But yes this set an arbitrary base resolution for the extent of the points, and build a pyramid up from there. It doesn't match how tyler works - Tyler has a fixed pyamid.

I think there will be a lot of use cases like this one (like just plotting big raster pyramids) so probably tyler needs to loosen up about where the zoom levels and tiles are meant to be.

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.