init code

Author: behinger

Project.toml

@@ -3,5 +3,19 @@ uuid = "a2bd6b23-9ae5-5b43-94df-c047a0516ae5"
 authors = ["Benedikt Ehinger <[email protected]>"]
 version = "0.1.0"
 
+[deps]
+Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
+Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
 [compat]
+Atomix = "1.0.1"
+Distances = "0.10.12"
+Distributions = "0.25.117"
+Interpolations = "0.15.1"
+KernelAbstractions = "0.9.33"
+Random = "1.11.0"
 julia = "1.10"

docs/Project.toml

@@ -1,11 +1,13 @@
-# Don't forget to run
-#
-#     pkg> dev ..
-#
 [deps]
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+DiscreteVoronoi = "01fd5c3e-b3b7-45fb-9d19-e82518019796"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 MakieStippling = "a2bd6b23-9ae5-5b43-94df-c047a0516ae5"
+MeshGrid = "ebf956a0-ef5e-43be-9fb1-27952996e635"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008"
 
 [compat]
 Documenter = "1"

docs/src/01-stippling.md

@@ -0,0 +1,70 @@
+```julia
+using MakieStippling
+using CUDA
+#using Interpolations
+using Random
+
+using CairoMakie
+using MeshGrid
+```
+
+create random data
+
+```julia
+n_sites = 5000
+n_conditions = 2
+sites_grouping = rand(1:n_conditions,n_sites)
+
+
+grid = Int.(zeros(500,500))
+grid_combined = CuArray(grid)
+sites = CartesianIndices(size(grid))[rand(1:length(grid),n_sites)]
+nothing #hide
+```
+
+create an arbitrary density function
+
+```julia
+times = collect(range(0,size(grid_combined,1),length=5))
+density_interpolators = Vector{Any}(undef,n_conditions)
+max_density = 1
+for k = 1:n_conditions
+    mean_data = rand(MersenneTwister(k),5).*size(grid,1)/4 .+ size(grid,1)/2
+    std_data = [20,10,5,15,30.]#abs.(rand(MersenneTwister(k*2),5)).*size(grid,1)
+
+    itp_mean = MakieStippling.Interpolations.interpolate((times,), mean_data, MakieStippling.Interpolations.Gridded(MakieStippling.Interpolations.Linear()))
+    itp_std = MakieStippling.Interpolations.interpolate((times,), std_data, MakieStippling.Interpolations.Gridded(MakieStippling.Interpolations.Linear()))
+
+    density_interpolators[k] = MakieStippling.DensityInterpolator(itp_mean,itp_std,max_density)
+
+end
+
+
+x,y = meshgrid(1:size(grid_combined,1),1:size(grid_combined,2))
+pts = Point2f.(x,y)
+
+fig = Figure()
+ax_dens = fig[1,1] = Axis(fig,title="underlying densities")
+hidespines!(ax_dens)
+hidedecorations!(ax_dens)
+for k = 1:n_conditions
+heatmap(fig[1,1][1,k],MakieStippling.calculate_density.(Ref(density_interpolators[k]),pts)')
+hidedecorations!(current_axis())
+end
+fig
+```
+
+```julia
+grid_single,sites_sets = MakieStippling.run_iterations!(grid_combined,sites,sites_grouping,density_interpolators;n_iter=3,threshold_low=10e-10,threshold_high=0.65*10-6)
+
+
+site_set_to_point2f(s::Set) = [Point2f(k[1],k[2]) for k in MakieStippling.linearized_sites.(collect(s),Ref(size(grid_combined)))]
+
+fig = Figure()
+ax = fig[1,1] = Axis(fig,aspect=1)
+scatter!.(ax,site_set_to_point2f.(sites_sets);markerspace = :data,)
+xlims!(ax,[0,500])
+ylims!(ax,[0,500])
+fig
+
+```

docs/src/02-comparing jfa implementations.md

@@ -0,0 +1,60 @@
+```julia
+using MakieStippling
+using CUDA
+using DiscreteVoronoi
+using StaticArrays
+using CairoMakie
+
+
+```
+
+```julia
+grid = Int.(zeros(30,30))
+grid_static = zeros(SVector{2,Int}, size(grid,1),size(grid,2))
+cu_grid = CuArray(grid)
+
+n_sites = 100
+sites = CartesianIndices(size(grid))[rand(1:length(grid),n_sites)]
+#sites_lin = MakieStippling.linearized_sites(sites, size(grid))
+
+sites_static = [SVector{2,Int}(x,y) for (x,y) in Tuple.(sites)]
+#sites = [CartesianIndex(1,1),CartesianIndex(10,10)]
+#grid[sites] .= 1:size(sites,1)
+#---
+MakieStippling.jfa_voronoi_gpu!(grid,sites)
+MakieStippling.jfa_voronoi_gpu!(cu_grid,sites)
+
+DiscreteVoronoi.jfa_voronoi!(grid_static,sites_static)
+
+# convert to linearized
+grid_static_lin =MakieStippling.linearized_sites(grid_static,size(grid_static))
+
+n_elems = *(size(grid)...)
+
+# DiscreteVoronoi vs. CPU
+ sum(grid_static_lin .== grid) / n_elems
+```
+
+```julia
+# CPU vs GPU
+sum(Matrix(cu_grid) .== grid) / n_elems
+```
+
+```julia
+# DiscreteVoronoi vs. GPU
+sum(Matrix(cu_grid) .== grid_static_lin)/n_elems
+
+```
+
+```julia
+f = Figure()
+
+heatmap(f[1,1],grid_static_lin,colormap=:tab10,axis=(;title="DiscreteVoronoi.jl"))
+heatmap(f[1,2],grid,colormap=:tab10,axis=(;title="cpu"))
+heatmap(f[1,3],Matrix(cu_grid),colormap=:tab10,axis=(;title="cu"))
+
+heatmap(f[2,1],grid.==grid_static_lin,axis=(;title="grid - DiscreteVoronoi.jl"))
+heatmap(f[2,2],Matrix(cu_grid).==grid_static_lin,axis=(;title="Cu - DiscreteVoronoi.jl"))
+heatmap(f[2,3],Matrix(cu_grid).==grid,axis=(;title="Cu - Grid"))
+
+f

src/MakieStippling.jl

@@ -1,11 +1,19 @@
 module MakieStippling
+using Distributions: value_support
 
-"""
-    hi = hello_world()
-A simple function to return "Hello, World!"
-"""
-function hello_world()
-    return "Hello, World!"
-end
+using Interpolations, Random
+using Distributions
+
+using Distances
+using KernelAbstractions
+
+using Atomix
+
+
+
+include("jfa_voronoi.jl")
+include("grid_features.jl")
+include("density.jl")
+include("stippling.jl")
 
 end

src/density.jl

@@ -0,0 +1,17 @@
+
+
+
+struct DensityInterpolator
+    mean_interpolator::Any
+    std_interpolator::Any
+    max_density::Any
+end
+
+function calculate_density(d::DensityInterpolator, point)
+    #@show point
+    mean = d.mean_interpolator(point[1])
+    std = d.std_interpolator(point[1])
+    pdf_value = pdf(Normal(mean, std), point[2])
+    normalized_value = (pdf_value) / (d.max_density)
+    return normalized_value
+end

src/grid_features.jl

@@ -0,0 +1,102 @@
+
+
+"""
+returns area and centroids for all elements in sites_set.
+"""
+function grid_features(A::AbstractMatrix{<:Integer}, sites_set)
+    device = get_backend(A)
+    # Define the number of bins for the histogram
+    num_bins = length(sites_set)
+
+    # Create arrays to store the histogram counts
+    hist_counts = similar(A, Int32, num_bins)
+    hist_counts .= 0
+
+    # Create arrays to store the sum of x-coordinates for each bin
+    hist_x_sum = similar(A, Int32, num_bins)
+    hist_x_sum .= 0
+    # Create arrays to store the sum of y-coordinates for each bin
+    hist_y_sum = similar(A, Int32, num_bins)
+    hist_y_sum .= 0
+
+    sites_array = similar(A, Int32, num_bins)
+    sites_array = typeof(sites_array)(Int32.(collect(sites_set)))
+    # Define the kernel
+    @kernel function histogram_kernel!(A, hist_counts, hist_x_sum, hist_y_sum, sites_set)
+        i = @index(Global, Linear)
+        #x = mod(i - 1, size(A, 2)) + 1
+        #y = cld(i, size(A, 2))
+        sz_grid = size(A)
+        coords = CartesianIndices(sz_grid)
+
+        x = coords[i][1]
+        y = coords[i][2]
+
+        bin_idx = findfirst(==(A[x, y]), sites_set)
+        if !isnothing(bin_idx)
+            Atomix.@atomic hist_counts[bin_idx] += 1
+            KernelAbstractions.@atomic hist_x_sum[bin_idx] += Int32(x)
+            KernelAbstractions.@atomic hist_y_sum[bin_idx] += Int32(y)
+            #=
+            =#
+        end
+    end
+    # @show "launch hist"
+    # Launch the kernel
+    event = histogram_kernel!(device, 128)(
+        A,
+        hist_counts,
+        hist_x_sum,
+        hist_y_sum,
+        sites_array;
+        ndrange = size(A),
+    )
+    KernelAbstractions.synchronize(device)
+
+
+    # Calculate the centroids
+    centroids_x = similar(A, Float32, num_bins)
+    centroids_y = similar(A, Float32, num_bins)
+
+    @kernel function centroid_kernel!(
+        hist_counts,
+        hist_x_sum,
+        hist_y_sum,
+        centroids_x,
+        centroids_y,
+    )
+        i = @index(Global, Linear)
+        if hist_counts[i] > 0
+            centroids_x[i] = hist_x_sum[i] / hist_counts[i]
+            centroids_y[i] = hist_y_sum[i] / hist_counts[i]
+        end
+    end
+    #    @show "launch centroid"
+    # Launch the centroid kernel
+    event = centroid_kernel!(device, 128)(
+        hist_counts,
+        hist_x_sum,
+        hist_y_sum,
+        centroids_x,
+        centroids_y;
+        ndrange = (num_bins,),
+    )
+    KernelAbstractions.synchronize(device)
+
+    return Vector(hist_counts) ./ *(size(A)...),
+    vec2tuple.(Vector(centroids_x), Vector(centroids_y))
+end
+vec2tuple(x, y) = (x, y)
+
+
+
+#=
+naive CPU version
+
+function grid_features(grid,id::Int)
+    ix = findall(x->x==id,grid)
+    centroid = Tuple(sum(ix))./length(ix)
+    area = length(ix)
+return centroid,area
+end
+=#

src/jfa_voronoi.jl

@@ -0,0 +1,60 @@
+
+function jfa_voronoi_gpu!(grid, sites::AbstractVector{<:CartesianIndex}; kwargs...)
+    sz_grid = size(grid)
+    jfa_voronoi_gpu!(grid, linearized_sites(sites, sz_grid); kwargs...)
+end
+
+function jfa_voronoi_gpu!(
+    grid,
+    sites::Vector{Int};
+    distance = euclidean,
+    do_1plusjfa = false,
+)
+    grid .= 0
+    grid[sites] = sites
+    k = max(size(grid)...)
+
+    if do_1plusjfa
+
+        jfa_kernel(get_backend(grid), 128)(grid, 1, ndrange = size(grid))
+        KernelAbstractions.synchronize(get_backend(grid))
+    end
+
+    while k > 1
+        k = k ÷ 2 + k % 2
+        jfa_kernel(get_backend(grid), 128)(grid, k, distance, ndrange = size(grid))
+        KernelAbstractions.synchronize(get_backend(grid))
+    end
+end
+
+
+
+@kernel function jfa_kernel(grid, k, distance)
+    ix_grid = @index(Global) # could be slightly optimized with Global, Cartesian
+
+    sz_grid = size(grid)
+    coords = CartesianIndices(sz_grid)
+
+    coords_grid = Tuple(coords[ix_grid])
+    (x, y) = coords_grid
+
+    for j in (-k, 0, k), i in (-k, 0, k)
+        checkbounds(Bool, grid, x + i, y + j) || continue
+        i == j == 0 && continue
+        # get candidate at k-distanced grid location
+        ix_siteq = grid[x+i, y+j]
+        ix_siteq != 0 || continue
+        # get current gridlocation value
+        ix_sitep = grid[x, y]
+        # if current location is not yet assigned, assign it with candidate
+        if ix_sitep == 0
+            grid[x, y] = ix_siteq
+        elseif distance(Tuple(coords[ix_sitep]), coords_grid) >
+               distance(Tuple(coords[ix_siteq]), coords_grid)
+            # if candidate site q is closer to what is currently saved, use that
+            grid[x, y] = ix_siteq
+        end
+    end
+
+
+end