Laplace.jl
edited
1function ∇²(A::Array{Float32,2})
2 """ Compute the Laplace Transform via Direct Convolution """
3 k_ℒ = Float32[[0.05 0.2 0.05]
4 [0.2 -1 0.2 ]
5 [0.05 0.2 0.05]]
6 A_prime = fill(0.0f0, size(A))
7
8 Threads.@threads for y = 1:size(A, 2)
9 Threads.@threads for x = 1:size(A,1)
10 # Bound A
11 by_A = max(1, y - 1) : min(size(A,2), y + 1)
12 bx_A = max(1, x - 1) : min(size(A,1), x + 1)
13
14 # Bound Laplace Kernel
15 by_ℒ = 1 + δ(y, 1) : size(k_ℒ, 2) - δ(y, size(A, 2))
16 bx_ℒ = 1 + δ(x, 1) : size(k_ℒ, 1) - δ(x, size(A, 1))
17
18 # Calculate subarray views
19 v_A = view(A, by_A, bx_A)
20 v_ℒ = view(k_ℒ, by_ℒ, bx_ℒ)
21
22 # Compute transform
23 A_prime[y, x] = sum(v_ℒ .* v_A)
24 end
25 end
26
27 return A_prime
28end