diff --git a/Manifest.toml b/Manifest.toml index a3c4013..344508d 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.8.0" manifest_format = "2.0" -project_hash = "a4d9eae2fa8abdf1a8b5713c28157ffc7cf118b0" +project_hash = "978ab5100e50fd0e20314354462cb53285a870ca" [[deps.AbstractFFTs]] deps = ["ChainRulesCore", "LinearAlgebra"] @@ -16,6 +16,12 @@ git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481" uuid = "6e696c72-6542-2067-7265-42206c756150" version = "1.1.4" +[[deps.Accessors]] +deps = ["Compat", "CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "LinearAlgebra", "MacroTools", "Requires", "Test"] +git-tree-sha1 = "8557017cfc7b58baea05a43ed35538857e6d35b4" +uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" +version = "0.1.19" + [[deps.Adapt]] deps = ["LinearAlgebra"] git-tree-sha1 = "195c5505521008abea5aee4f96930717958eac6f" @@ -282,12 +288,23 @@ git-tree-sha1 = "d7f719f3c186e4a821fda1ab38e9cfeaf3b17b71" uuid = "6cd66ae4-5932-4b96-926d-e73e578e42cc" version = "0.1.9" +[[deps.FoldsThreads]] +deps = ["Accessors", "FunctionWrappers", "InitialValues", "SplittablesBase", "Transducers"] +git-tree-sha1 = "eb8e1989b9028f7e0985b4268dabe94682249025" +uuid = "9c68100b-dfe1-47cf-94c8-95104e173443" +version = "0.1.1" + [[deps.ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"] git-tree-sha1 = "187198a4ed8ccd7b5d99c41b69c679269ea2b2d4" uuid = "f6369f11-7733-5829-9624-2563aa707210" version = "0.10.32" +[[deps.FunctionWrappers]] +git-tree-sha1 = "241552bc2209f0fa068b6415b1942cc0aa486bcc" +uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" +version = "1.1.2" + [[deps.Future]] deps = ["Random"] uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" diff --git a/Project.toml b/Project.toml index 2dfc5e4..78c2245 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.1.0" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" FLoops = "cc61a311-1640-44b5-9fba-1b764f453329" FoldsCUDA = "6cd66ae4-5932-4b96-926d-e73e578e42cc" +FoldsThreads = "9c68100b-dfe1-47cf-94c8-95104e173443" ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" diff --git a/src/squared_euclidean.jl b/src/squared_euclidean.jl index ec1a801..acbc373 100644 --- a/src/squared_euclidean.jl +++ b/src/squared_euclidean.jl @@ -15,6 +15,7 @@ begin using DistanceTransforms using FLoops using CUDA + using FoldsThreads end # ╔═╡ 87dd183b-a462-4744-a800-3936791aea43 @@ -201,7 +202,7 @@ elements. Multi-Threaded version of `transform!(..., tfm::SquaredEuclidean)` """ # ╔═╡ 4f4bd88d-9941-486c-b562-ec3fb9cadba7 -function transform!(img::AbstractMatrix, tfm::SquaredEuclidean, nthreads; output=zeros(size(img)), v=ones(Int32, size(img)), z=ones(size(img) .+ 1)) +function transform!(img::AbstractMatrix, tfm::SquaredEuclidean, nthreads::Number; output=zeros(size(img)), v=ones(Int32, size(img)), z=ones(size(img) .+ 1)) Threads.@threads for i in axes(img, 1) @views transform(img[i, :], tfm; output=output[i,:], v=fill!(v[i,:], 1), z=fill!(z[i,:], 1)) end @@ -222,7 +223,7 @@ Applies a squared euclidean distance transform to an input image. Returns an arr """ # ╔═╡ e570b1d5-b2ff-48c0-8562-9cf2d350f9e4 -function transform!(vol::AbstractArray, tfm::SquaredEuclidean, nthreads; output=zeros(size(vol)), v=ones(Int32, size(vol)), z=ones(size(vol) .+ 1)) +function transform!(vol::AbstractArray, tfm::SquaredEuclidean, nthreads::Number; output=zeros(size(vol)), v=ones(Int32, size(vol)), z=ones(size(vol) .+ 1)) Threads.@threads for k in axes(vol, 3) @views transform!(vol[:, :, k], tfm; output=output[:, :, k], v=fill!(v[:, :, k], 1), z=fill!(z[:, :, k], 1)) end @@ -234,17 +235,6 @@ function transform!(vol::AbstractArray, tfm::SquaredEuclidean, nthreads; output= return output end -# ╔═╡ e559548d-1234-489a-b59e-db80aaf81798 -# function transform!(vol::AbstractArray, tfm::SquaredEuclidean, nthreads; output=zeros(size(vol)), v=ones(Int32, size(vol)), z=ones(size(vol) .+ 1)) -# Threads.@threads for k in axes(vol, 3) -# @views transform!(vol[:, :, k], tfm; output=output[:, :, k], v=fill!(v[:, :, k], 1), z=fill!(z[:, :, k], 1)) -# end -# Threads.@threads for (i, j) in (axes(vol, 2), axes(vol, 1)) -# @views transform(output[i, j, :], tfm; output=output[i, j, :], v=fill!(v[i, j, :], 1), z=fill!(z[i, j, :], 1)) -# end -# return output -# end - # ╔═╡ a5a497b8-d1c6-4f00-8ab1-75dc9571cc0a md""" ## GPU @@ -284,6 +274,45 @@ function transform!(vol::CuArray{T, 3}, tfm::SquaredEuclidean; output=CUDA.zeros return output end +# ╔═╡ 3bb1c028-4f7a-414a-86bd-a4522af64957 +md""" +## Various Multi-Threading +""" + +# ╔═╡ d470c994-4de2-4d84-950f-4508b11827a2 +md""" +### 2D! +""" + +# ╔═╡ 0726f423-2044-4d78-8eca-9433e9f1dc95 +function transform!(img::AbstractMatrix, tfm::SquaredEuclidean, ex; output=CUDA.zeros(size(img)), v=CUDA.ones(size(img)), z=CUDA.ones(size(img) .+ 1)) where T + @floop ex for i in axes(img, 1) + @views transform(img[i, :], tfm; output=output[i,:], v=v[i,:], z=z[i,:]) + end + @floop ex for j in axes(img, 2) + @views transform(output[:, j], tfm; output=output[:,j], v=fill!(v[:,j], 1), z=fill!(z[:,j], 1)) + end + return output +end + +# ╔═╡ f2ceeaa5-988d-4adb-9ac2-67c0ded6198d +md""" +### 3D! +""" + +# ╔═╡ 07e9ebfa-a287-45d9-a38a-aab7690992f9 +function transform!(vol::AbstractArray, tfm::SquaredEuclidean, ex; output=zeros(size(vol)), v=ones(size(vol)), z=ones(size(vol) .+ 1)) where T + @floop ex for k in axes(vol, 3) + @views transform!(vol[:, :, k], tfm; output=output[:, :, k], v=v[:, :, k], z=z[:, :, k]) + end + @floop ex for i in axes(vol, 1) + for j in axes(vol, 2) + @views transform(output[i, j, :], tfm; output=output[i, j, :], v=fill!(v[i, j, :], 1), z=fill!(z[i, j, :], 1)) + end + end + return output +end + # ╔═╡ Cell order: # ╠═c7ef37d8-2330-11ed-006d-a16889a98cd1 # ╠═87dd183b-a462-4744-a800-3936791aea43 @@ -306,9 +335,13 @@ end # ╠═4f4bd88d-9941-486c-b562-ec3fb9cadba7 # ╟─1df8eb42-e73c-4680-a8ce-c9d648ee8fe4 # ╠═e570b1d5-b2ff-48c0-8562-9cf2d350f9e4 -# ╠═e559548d-1234-489a-b59e-db80aaf81798 # ╟─a5a497b8-d1c6-4f00-8ab1-75dc9571cc0a # ╟─2ba78c10-a3c2-460c-b30e-d86fc04c581f # ╠═b8a0a7a0-887f-4c87-80ad-41a28aa8bf1c # ╟─c5d95fe5-fb5f-441d-9392-1ac9eec06924 # ╠═4333808c-d4b5-479c-b659-78fc0c17bf51 +# ╟─3bb1c028-4f7a-414a-86bd-a4522af64957 +# ╟─d470c994-4de2-4d84-950f-4508b11827a2 +# ╠═0726f423-2044-4d78-8eca-9433e9f1dc95 +# ╟─f2ceeaa5-988d-4adb-9ac2-67c0ded6198d +# ╠═07e9ebfa-a287-45d9-a38a-aab7690992f9 diff --git a/test/squared_euclidean.jl b/test/squared_euclidean.jl index 8a47828..14acf38 100644 --- a/test/squared_euclidean.jl +++ b/test/squared_euclidean.jl @@ -14,6 +14,7 @@ begin using Test using CUDA using DistanceTransforms + using FoldsThreads end # ╔═╡ 373a8802-bbac-4ab5-abe4-420ccbb61ea0 @@ -467,7 +468,7 @@ if CUDA.has_cuda_gpu() @test test == CuArray(answer) end else - @warn "CUDA" + @warn "CUDA unavailable, not testing GPU support" end; # ╔═╡ 8998856b-dc41-4774-8e04-972caf278e10 @@ -487,6 +488,279 @@ else end; +# ╔═╡ 3e96bc74-9519-4547-b5d0-62e767a28b94 +md""" +## Various Multi-Threading +""" + +# ╔═╡ 7db96dd6-cddd-4cb6-8851-b682b64b6fb2 +md""" +### 2D! +""" + +# ╔═╡ e24c93b2-795a-45ab-a40d-215cbee22660 +@testset "squared euclidean 2D FoldsThreads " begin + img = [ + 0 1 1 1 0 0 0 1 1 + 1 1 1 1 1 0 0 0 1 + 1 0 0 0 1 0 0 1 1 + 1 0 0 0 1 0 1 1 0 + 1 0 0 0 1 1 0 1 0 + 1 1 1 1 1 0 0 1 0 + 0 1 1 1 0 0 0 0 1 + ] + output, v, z = zeros(size(img)), ones(Int32, size(img)), ones(size(img) .+ 1) + tfm = SquaredEuclidean() + ex = DepthFirstEx() + test = transform!(boolean_indicator(img), tfm, ex; output=output, v=v, z=z) + answer = [ + 1.0 0.0 0.0 0.0 1.0 2.0 1.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 1.0 2.0 1.0 0.0 + 0.0 1.0 1.0 1.0 0.0 1.0 1.0 0.0 0.0 + 0.0 1.0 4.0 1.0 0.0 1.0 0.0 0.0 1.0 + 0.0 1.0 1.0 1.0 0.0 0.0 1.0 0.0 1.0 + 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 1.0 + 1.0 0.0 0.0 0.0 1.0 2.0 2.0 1.0 0.0 + ] + @test test == answer +end; + +# ╔═╡ f44983e3-af78-4dc6-8447-2888806af34d +@testset "squared euclidean 2D FoldsThreads " begin + img = rand([0, 1], 10, 10) + output, v, z = zeros(size(img)), ones(Int32, size(img)), ones(size(img) .+ 1) + tfm = SquaredEuclidean() + ex = DepthFirstEx() + test = transform!(boolean_indicator(img), tfm, ex; output=output, v=v, z=z) + answer = transform(boolean_indicator(img), tfm) + @test test == answer +end; + +# ╔═╡ 7cdb1f83-1323-4237-bf3c-eac5b849a514 +@testset "squared euclidean 2D FoldsThreads " begin + img = [ + 0 1 1 1 0 0 0 1 1 + 1 1 1 1 1 0 0 0 1 + 1 0 0 0 1 0 0 1 1 + 1 0 0 0 1 0 1 1 0 + 1 0 0 0 1 1 0 1 0 + 1 1 1 1 1 0 0 1 0 + 0 1 1 1 0 0 0 0 1 + ] + output, v, z = zeros(size(img)), ones(Int32, size(img)), ones(size(img) .+ 1) + tfm = SquaredEuclidean() + ex = NonThreadedEx() + test = transform!(boolean_indicator(img), tfm, ex; output=output, v=v, z=z) + answer = [ + 1.0 0.0 0.0 0.0 1.0 2.0 1.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 1.0 2.0 1.0 0.0 + 0.0 1.0 1.0 1.0 0.0 1.0 1.0 0.0 0.0 + 0.0 1.0 4.0 1.0 0.0 1.0 0.0 0.0 1.0 + 0.0 1.0 1.0 1.0 0.0 0.0 1.0 0.0 1.0 + 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 1.0 + 1.0 0.0 0.0 0.0 1.0 2.0 2.0 1.0 0.0 + ] + @test test == answer +end; + +# ╔═╡ 17ee7bbf-8220-4549-8458-1723e5a3e30d +@testset "squared euclidean 2D FoldsThreads " begin + img = rand([0, 1], 10, 10) + output, v, z = zeros(size(img)), ones(Int32, size(img)), ones(size(img) .+ 1) + tfm = SquaredEuclidean() + ex = NonThreadedEx() + test = transform!(boolean_indicator(img), tfm, ex; output=output, v=v, z=z) + answer = transform(boolean_indicator(img), tfm) + @test test == answer +end; + +# ╔═╡ 74e9c7b1-ad17-4759-895d-550133b09788 +@testset "squared euclidean 2D FoldsThreads " begin + img = [ + 0 1 1 1 0 0 0 1 1 + 1 1 1 1 1 0 0 0 1 + 1 0 0 0 1 0 0 1 1 + 1 0 0 0 1 0 1 1 0 + 1 0 0 0 1 1 0 1 0 + 1 1 1 1 1 0 0 1 0 + 0 1 1 1 0 0 0 0 1 + ] + output, v, z = zeros(size(img)), ones(Int32, size(img)), ones(size(img) .+ 1) + tfm = SquaredEuclidean() + ex = WorkStealingEx() + test = transform!(boolean_indicator(img), tfm, ex; output=output, v=v, z=z) + answer = [ + 1.0 0.0 0.0 0.0 1.0 2.0 1.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 1.0 2.0 1.0 0.0 + 0.0 1.0 1.0 1.0 0.0 1.0 1.0 0.0 0.0 + 0.0 1.0 4.0 1.0 0.0 1.0 0.0 0.0 1.0 + 0.0 1.0 1.0 1.0 0.0 0.0 1.0 0.0 1.0 + 0.0 0.0 0.0 0.0 0.0 1.0 1.0 0.0 1.0 + 1.0 0.0 0.0 0.0 1.0 2.0 2.0 1.0 0.0 + ] + @test test == answer +end; + +# ╔═╡ 9f78df21-d298-495f-a656-5325807cb109 +@testset "squared euclidean 2D FoldsThreads" begin + img = rand([0, 1], 10, 10) + output, v, z = zeros(size(img)), ones(Int32, size(img)), ones(size(img) .+ 1) + tfm = SquaredEuclidean() + ex = WorkStealingEx() + test = transform!(boolean_indicator(img), tfm, ex; output=output, v=v, z=z) + answer = transform(boolean_indicator(img), tfm) + @test test == answer +end; + +# ╔═╡ ed548af6-f1b9-4302-b787-b8ad67d4a5b4 +md""" +### 3D! +""" + +# ╔═╡ e19e8cc6-1f9c-4fa8-9783-e148b6729bb7 +@testset "squared euclidean 3D FoldsThreads" begin + img = [ + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 1 1 1 0 0 0 0 0 + 0 0 1 0 0 1 0 0 0 0 0 + 0 0 1 0 0 1 1 1 0 0 0 + 0 0 1 0 0 0 0 1 0 0 0 + 0 0 1 0 0 0 0 1 0 0 0 + 0 0 0 1 1 1 1 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + ] + img_inv = @. ifelse(img == 0, 1, 0) + vol = cat(img, img_inv, dims=3) + container2 = [] + for i in 1:10 + push!(container2, vol) + end + vol_inv = cat(container2..., dims=3) + output, v, z = zeros(size(vol_inv)), ones(Int32, size(vol_inv)), ones(size(vol_inv) .+ 1) + tfm = SquaredEuclidean() + ex = DepthFirstEx() + test = transform!(boolean_indicator(vol_inv), tfm, ex; output=output, v=v, z=z) + a1 = img_inv + a2 = img + ans = cat(a1, a2, dims=3) + container_a = [] + for i in 1:10 + push!(container_a, ans) + end + answer = cat(container_a..., dims=3) + @test test == answer +end; + +# ╔═╡ 8c570c8b-ffc4-4b2b-845a-a504b179af4e +@testset "squared euclidean 3D FoldsThreads" begin + img = rand([0, 1], 10, 10, 10) + output, v, z = zeros(size(img)), ones(Int32, size(img)), ones(size(img) .+ 1) + tfm = SquaredEuclidean() + ex = WorkStealingEx() + test = transform!(boolean_indicator(img), tfm, ex; output=output, v=v, z=z) + answer = transform(boolean_indicator(img), tfm) + @test test == answer +end; + +# ╔═╡ a3c842dc-f392-4505-9cf2-3ab54cd80955 +@testset "squared euclidean 3D FoldsThreads" begin + img = [ + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 1 1 1 0 0 0 0 0 + 0 0 1 0 0 1 0 0 0 0 0 + 0 0 1 0 0 1 1 1 0 0 0 + 0 0 1 0 0 0 0 1 0 0 0 + 0 0 1 0 0 0 0 1 0 0 0 + 0 0 0 1 1 1 1 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + ] + img_inv = @. ifelse(img == 0, 1, 0) + vol = cat(img, img_inv, dims=3) + container2 = [] + for i in 1:10 + push!(container2, vol) + end + vol_inv = cat(container2..., dims=3) + output, v, z = zeros(size(vol_inv)), ones(Int32, size(vol_inv)), ones(size(vol_inv) .+ 1) + tfm = SquaredEuclidean() + ex = NonThreadedEx() + test = transform!(boolean_indicator(vol_inv), tfm, ex; output=output, v=v, z=z) + a1 = img_inv + a2 = img + ans = cat(a1, a2, dims=3) + container_a = [] + for i in 1:10 + push!(container_a, ans) + end + answer = cat(container_a..., dims=3) + @test test == answer +end; + +# ╔═╡ 9850ad92-53cb-431d-92c0-04d61c4db391 +@testset "squared euclidean 3D FoldsThreads" begin + img = rand([0, 1], 10, 10, 10) + output, v, z = zeros(size(img)), ones(Int32, size(img)), ones(size(img) .+ 1) + tfm = SquaredEuclidean() + ex = NonThreadedEx() + test = transform!(boolean_indicator(img), tfm, ex; output=output, v=v, z=z) + answer = transform(boolean_indicator(img), tfm) + @test test == answer +end; + +# ╔═╡ 1d7f7eb8-2b7d-4d5c-aa09-bdb324047443 +@testset "squared euclidean 3D FoldsThreads" begin + img = [ + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 1 1 1 0 0 0 0 0 + 0 0 1 0 0 1 0 0 0 0 0 + 0 0 1 0 0 1 1 1 0 0 0 + 0 0 1 0 0 0 0 1 0 0 0 + 0 0 1 0 0 0 0 1 0 0 0 + 0 0 0 1 1 1 1 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 + ] + img_inv = @. ifelse(img == 0, 1, 0) + vol = cat(img, img_inv, dims=3) + container2 = [] + for i in 1:10 + push!(container2, vol) + end + vol_inv = cat(container2..., dims=3) + output, v, z = zeros(size(vol_inv)), ones(Int32, size(vol_inv)), ones(size(vol_inv) .+ 1) + tfm = SquaredEuclidean() + ex = WorkStealingEx() + test = transform!(boolean_indicator(vol_inv), tfm, ex; output=output, v=v, z=z) + a1 = img_inv + a2 = img + ans = cat(a1, a2, dims=3) + container_a = [] + for i in 1:10 + push!(container_a, ans) + end + answer = cat(container_a..., dims=3) + @test test == answer +end; + +# ╔═╡ fe5fb42c-9c56-4c27-a2c5-70b8e761d9b8 +@testset "squared euclidean 3D FoldsThreads" begin + img = rand([0, 1], 10, 10, 10) + output, v, z = zeros(size(img)), ones(Int32, size(img)), ones(size(img) .+ 1) + tfm = SquaredEuclidean() + ex = WorkStealingEx() + test = transform!(boolean_indicator(img), tfm, ex; output=output, v=v, z=z) + answer = transform(boolean_indicator(img), tfm) + @test test == answer +end; + # ╔═╡ Cell order: # ╠═cd8bf944-2329-11ed-208f-1b2e91673a5e # ╠═373a8802-bbac-4ab5-abe4-420ccbb61ea0 @@ -522,3 +796,18 @@ end; # ╟─b5ba030e-9d5e-4ce5-b138-f3be7fea0e23 # ╠═223acdc8-b5de-4eb1-9c80-19c3281e0dfd # ╠═8998856b-dc41-4774-8e04-972caf278e10 +# ╟─3e96bc74-9519-4547-b5d0-62e767a28b94 +# ╟─7db96dd6-cddd-4cb6-8851-b682b64b6fb2 +# ╠═e24c93b2-795a-45ab-a40d-215cbee22660 +# ╠═f44983e3-af78-4dc6-8447-2888806af34d +# ╠═7cdb1f83-1323-4237-bf3c-eac5b849a514 +# ╠═17ee7bbf-8220-4549-8458-1723e5a3e30d +# ╠═74e9c7b1-ad17-4759-895d-550133b09788 +# ╠═9f78df21-d298-495f-a656-5325807cb109 +# ╟─ed548af6-f1b9-4302-b787-b8ad67d4a5b4 +# ╠═e19e8cc6-1f9c-4fa8-9783-e148b6729bb7 +# ╠═8c570c8b-ffc4-4b2b-845a-a504b179af4e +# ╠═a3c842dc-f392-4505-9cf2-3ab54cd80955 +# ╠═9850ad92-53cb-431d-92c0-04d61c4db391 +# ╠═1d7f7eb8-2b7d-4d5c-aa09-bdb324047443 +# ╠═fe5fb42c-9c56-4c27-a2c5-70b8e761d9b8