Skip to content

Instantly share code, notes, and snippets.

@rafaqz
Last active April 30, 2023 03:20
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save rafaqz/a1706c16418bf875cc3aafc631dda49a to your computer and use it in GitHub Desktop.
"""
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
Copy link

asinghvi17 commented Apr 29, 2023

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
Copy link
Author

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
Copy link

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment