Julia Implementation of Integral Arrays
IntegralArrays are useful for summing arrays over rectangular regions. Once created, the sum of an arbitrarily-large rectangular region
can be computed in O(2^d) adds, where d is the dimensionality.
Demo:
julia> using IntegralArrays, IntervalSets
julia> A = [1 2 3 4 5;
10 20 30 40 50;
100 200 300 400 500]
3×5 Matrix{Int64}:
1 2 3 4 5
10 20 30 40 50
100 200 300 400 500
julia> Ai = IntegralArray(A)
3×5 IntegralArray{Int64, 2, Matrix{Int64}}:
1 3 6 10 15
11 33 66 110 165
111 333 666 1110 1665
julia> sum(A[1:2, 1:2]) == Ai[2, 2]
true
julia> sum(A[1:3, 2:4]) == Ai[1..3, 2..4]
trueWhen one needs to compute the sum/average of all blocks extracted from an image, pre-building the integral array usually provides a more efficient computation.
using BenchmarkTools, IntegralArrays, IntervalSets
# simplified 3x3 mean filter; only for demo purpose
function mean_filter_naive!(out, X)
Δ = CartesianIndex(1, 1)
for i in CartesianIndex(2, 2):CartesianIndex(size(X).-1)
block = @view X[i-Δ: i+Δ]
out[i] = mean(block)
end
return out
end
function mean_filter_integral!(out, X)
iX = IntegralArray(X)
for i in CartesianIndex(2, 2):CartesianIndex(size(X).-1)
x, y = i.I
out[i] = iX[x±1, y±1]/9
end
return out
end
X = Float32.(rand(1:5, 64, 64));
m1 = copy(X);
m2 = copy(X);
@btime mean_filter_naive!($m1, $X); # 65.078 μs (0 allocations: 0 bytes)
@btime mean_filter_integral!($m2, $X); # 12.161 μs (4 allocations: 16.17 KiB)
m1 == m2 # true