Instantly share code, notes, and snippets.

Last active April 30, 2023 03:20
Star You must be signed in to star a gist
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 """ 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)

### asinghvi17 commented Apr 29, 2023 • edited

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.

### rafaqz commented Apr 29, 2023

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.

### asinghvi17 commented Apr 30, 2023

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.