From a9bd4e40806b9eb8fa19523a226a9c7fd0e1f0df Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 18 Jan 2023 14:41:00 -0800 Subject: [PATCH 1/5] modify findfirst to return index rather than shape --- src/util/kdtree.jl | 32 +++++++++++++++++++++++++------- test/kdtree.jl | 2 +- 2 files changed, 26 insertions(+), 8 deletions(-) diff --git a/src/util/kdtree.jl b/src/util/kdtree.jl index 37f0669..7f0536e 100644 --- a/src/util/kdtree.jl +++ b/src/util/kdtree.jl @@ -16,14 +16,16 @@ not shapes of nonzero size.) """ mutable struct KDTree{K,S<:Shape{K}} s::Vector{S} + s_index::Vector{Int} ix::Int x::Float left::KDTree{K,S} # shapes ≤ x in coordinate ix right::KDTree{K,S} # shapes > x in coordinate ix - KDTree{K,S}(s::AbsVec{S}) where {K,S<:Shape{K}} = new(s, 0) + KDTree{K,S}(s::AbsVec{S}) where {K,S<:Shape{K}} = new(s, collect(1:size(s)[1]), 0) + KDTree{K,S}(s::AbsVec{S},s_index::Vector{<:Int}) where {K,S<:Shape{K}} = new(s, s_index, 0) function KDTree{K,S}(ix::Integer, x::Real, left::KDTree{K,S}, right::KDTree{K,S}) where {K,S<:Shape{K}} 1 ≤ ix ≤ K || throw(BoundsError()) - new(S[], ix, x, left, right) + new(S[], Int[], ix, x, left, right) end end @@ -38,8 +40,15 @@ Construct a K-D tree (`KDTree`) representation of a list of When searching the tree, shapes that appear earlier in `s` take precedence over shapes that appear later. """ + function KDTree(s::AbsVec{S}) where {K,S<:Shape{K}} - (length(s) ≤ 4 || K == 0) && return KDTree{K,S}(s) + # If no list of indicies is provided, simply enumerate by the number of + # shapes in `s`. + return KDTree(s,collect(1:size(s)[1])) +end + +function KDTree(s::AbsVec{S}, s_index::Vector{<:Int}) where {K,S<:Shape{K}} + (length(s) ≤ 4 || K == 0) && return KDTree{K,S}(s, s_index) # figure out the best dimension ix to divide over, # the dividing plane x, and the number (nl,nr) of @@ -61,22 +70,26 @@ function KDTree(s::AbsVec{S}) where {K,S<:Shape{K}} end # don't bother subdividing if it doesn't reduce the # of shapes much - 4*max(nl,nr) > 3*length(s) && return KDTree{K,S}(s) + 4*max(nl,nr) > 3*length(s) && return KDTree{K,S}(s,s_index) # create the arrays of shapes in each subtree sl = Vector{S}(undef, nl) + sl_idx = Vector{Int}(undef, nl) sr = Vector{S}(undef, nr) + sr_idx = Vector{Int}(undef, nr) il = ir = 0 for k in eachindex(s) if b[k][1][ix] ≤ x sl[il += 1] = s[k] + sl_idx[il] = s_index[k] end if b[k][2][ix] > x sr[ir += 1] = s[k] + sr_idx[ir] = s_index[k] end end - return KDTree{K,S}(ix, x, KDTree(sl), KDTree(sr)) + return KDTree{K,S}(ix, x, KDTree(sl,sl_idx), KDTree(sr,sr_idx)) end depth(kd::KDTree) = kd.ix == 0 ? 0 : max(depth(kd.left), depth(kd.right)) + 1 @@ -104,7 +117,7 @@ function Base.findfirst(p::SVec{N}, s::Vector{S}) where {N,S<:Shape{N}} for i in eachindex(s) b = bounds(s[i]) if all(b[1] .< p .< b[2]) && p ∈ s[i] # check if p is within bounding box is faster - return s[i] + return i end end return nothing @@ -118,7 +131,12 @@ function Base.findfirst(p::SVec{N}, kd::KDTree{N}) where {N} return findfirst(p, kd.right) end else - return findfirst(p, kd.s) + idx = findfirst(p, kd.s) + if isnothing(idx) + return idx + else + return kd.s_index[idx] + end end end diff --git a/test/kdtree.jl b/test/kdtree.jl index eab9209..9013e3c 100644 --- a/test/kdtree.jl +++ b/test/kdtree.jl @@ -7,7 +7,7 @@ s = Shape2[Ball([i,0], 1) for i in 0:20] kd = KDTree(s) @test GeometryPrimitives.depth(kd) == 3 - @test findfirst([10.1,0], kd).c[1] == 10 + @test s[findfirst([10.1,0], kd)].c[1] == 10 @test findfirst([10.1,1], kd) == nothing @test checktree(kd, s) From e8bce9ef7f73ef63a306adb9640c3ef0a8650149 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Thu, 19 Jan 2023 11:09:38 -0800 Subject: [PATCH 2/5] use eachindex instead of 1:size --- src/util/kdtree.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/util/kdtree.jl b/src/util/kdtree.jl index 7f0536e..4a27328 100644 --- a/src/util/kdtree.jl +++ b/src/util/kdtree.jl @@ -21,7 +21,7 @@ mutable struct KDTree{K,S<:Shape{K}} x::Float left::KDTree{K,S} # shapes ≤ x in coordinate ix right::KDTree{K,S} # shapes > x in coordinate ix - KDTree{K,S}(s::AbsVec{S}) where {K,S<:Shape{K}} = new(s, collect(1:size(s)[1]), 0) + KDTree{K,S}(s::AbsVec{S}) where {K,S<:Shape{K}} = new(s, collect(eachindex(s)), 0) KDTree{K,S}(s::AbsVec{S},s_index::Vector{<:Int}) where {K,S<:Shape{K}} = new(s, s_index, 0) function KDTree{K,S}(ix::Integer, x::Real, left::KDTree{K,S}, right::KDTree{K,S}) where {K,S<:Shape{K}} 1 ≤ ix ≤ K || throw(BoundsError()) @@ -44,7 +44,7 @@ take precedence over shapes that appear later. function KDTree(s::AbsVec{S}) where {K,S<:Shape{K}} # If no list of indicies is provided, simply enumerate by the number of # shapes in `s`. - return KDTree(s,collect(1:size(s)[1])) + return KDTree(s,collect(eachindex(s))) end function KDTree(s::AbsVec{S}, s_index::Vector{<:Int}) where {K,S<:Shape{K}} From 15e38b49a81aa3a98fa0cb9bb10935899b2deffe Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Fri, 20 Jan 2023 10:25:30 -0800 Subject: [PATCH 3/5] Use `Integer` not `Int` Co-authored-by: Steven G. Johnson --- src/util/kdtree.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/util/kdtree.jl b/src/util/kdtree.jl index 4a27328..b62c995 100644 --- a/src/util/kdtree.jl +++ b/src/util/kdtree.jl @@ -47,7 +47,7 @@ function KDTree(s::AbsVec{S}) where {K,S<:Shape{K}} return KDTree(s,collect(eachindex(s))) end -function KDTree(s::AbsVec{S}, s_index::Vector{<:Int}) where {K,S<:Shape{K}} +function KDTree(s::AbsVec{S}, s_index::AbstractVector{<:Integer}) where {K,S<:Shape{K}} (length(s) ≤ 4 || K == 0) && return KDTree{K,S}(s, s_index) # figure out the best dimension ix to divide over, From c13df0f7d4ea8becb0692bcdbf2ceab1063a77f3 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 26 Apr 2023 11:20:02 -0700 Subject: [PATCH 4/5] add simple smoothing example --- examples/simple_smoothing.jl | 61 ++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 examples/simple_smoothing.jl diff --git a/examples/simple_smoothing.jl b/examples/simple_smoothing.jl new file mode 100644 index 0000000..f9affed --- /dev/null +++ b/examples/simple_smoothing.jl @@ -0,0 +1,61 @@ +"""simple_smoothing.jl + +This example shows how one can use the underlying `surfpt_nearby()` and +`volfrac()` routines to perform simple subpixel smoothing (in this case, a +linear interpolation between a cylinder and a background void). + +Since we only assume a single geometric shape lies in the domain, we already +know what the foreground and background shapes are, and thus don't need to +search the tree at each pixel. If we did have a more complicated geometry +stackup, we would need to carefully determine what the foreground and background +materials are within each voxel (there could be more than 2 materials, in which +case there's no clear convention). +""" +using GeometryPrimitives +using CairoMakie +using StaticArrays + +function populate_matrix() + # set up geometry + geometry = [Cylinder([2.0, -2.0, 0.0], 1.0, 100.0)] + + # arbitrary material parameters + low = 1.0 + high = 3.0 + + # set up domain + sx, sy = 10.0, 10.0 + resolution = 5 # voxels/unit + + # build matrix + Nx = Int(sx * resolution) + Ny = Int(sy * resolution) + vxl_size = SVector(1.0 / resolution, 1.0 / resolution, 0.0) + M = zeros(Nx, Ny) + + # loop over grid points and fill the matrix + x = range(-sx / 2, stop = sx / 2, length = Nx) + y = range(-sy / 2, stop = sy / 2, length = Ny) + for ix = 1:Nx, iy = 1:Ny + p = SVector(x[ix], y[iy], 0.0) # current voxel center + vxl = (p - vxl_size / 2.0, p + vxl_size / 2.0) # current voxel on the grid + + # We only have one shape, so no need to "search" for any other shapes. + # But normally, one would want to search for a "foreground" shape and a + # "background" shape (within every pixel) and do averaging on those. + nearby = surfpt_nearby(p, geometry[1]) + # compute the fill fraction + fill_fraction = volfrac(vxl, nearby[2], nearby[1]) + # linearly interpolate based on the fill fraction + M[ix, iy] = low_index + (fill_fraction) * (high_index - low_index) + end + + return M +end + +function run_and_plot() + M = populate_matrix() + display(heatmap(M)) +end + +run_and_plot() From 1ca3d74d3c94d303754f42521e48a2159df06123 Mon Sep 17 00:00:00 2001 From: Alec Hammond Date: Wed, 26 Apr 2023 11:32:45 -0700 Subject: [PATCH 5/5] bug fixes --- examples/simple_smoothing.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/simple_smoothing.jl b/examples/simple_smoothing.jl index f9affed..cffd50c 100644 --- a/examples/simple_smoothing.jl +++ b/examples/simple_smoothing.jl @@ -11,13 +11,14 @@ stackup, we would need to carefully determine what the foreground and background materials are within each voxel (there could be more than 2 materials, in which case there's no clear convention). """ + using GeometryPrimitives -using CairoMakie +using Makie using StaticArrays function populate_matrix() # set up geometry - geometry = [Cylinder([2.0, -2.0, 0.0], 1.0, 100.0)] + geometry = [Cylinder([2.0, -3.0, 0.0], 1.0, 100.0)] # arbitrary material parameters low = 1.0 @@ -47,7 +48,7 @@ function populate_matrix() # compute the fill fraction fill_fraction = volfrac(vxl, nearby[2], nearby[1]) # linearly interpolate based on the fill fraction - M[ix, iy] = low_index + (fill_fraction) * (high_index - low_index) + M[ix, iy] = low + (fill_fraction) * (high - low) end return M