adds integral arrays. Closes JuliaImages/ImagesFiltering.jl#4 #1
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1 +1,3 @@ | ||
| julia 0.5- | ||
| Images | ||
| IntervalSets |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,5 +1,73 @@ | ||
| module IntegralArrays | ||
|
|
||
| # package code goes here | ||
| using Images, IntervalSets | ||
|
|
||
| import Base: size, getindex, linearindexing | ||
|
|
||
| export IntegralArray | ||
|
|
||
| """ | ||
| ``` | ||
| integral_array = IntegralArray(A) | ||
| ``` | ||
|
|
||
| Returns the integral array of an array. The integral array is calculated by assigning | ||
| to each cell the sum of all cells above it and to its left, i.e. the rectangle from | ||
| (1, 1) to the pixel. An integral array is a data structure which helps in efficient | ||
| calculation of sum of pixels in a rectangular subset of an array. | ||
|
|
||
| ``` | ||
| sum = integral_array[ytop..ybot, xtop..xbot] | ||
| sum = integral_array[y ± σy, x ± σx] | ||
| sum = integral_array[CartesianIndex(y, x) ± σ] | ||
| ``` | ||
|
|
||
| The sum of a window in an array can be directly calculated using four array references of | ||
| the integral array, irrespective of the size of the window, given the `yrange` and `xrange` | ||
| of the window. Given an integral array - | ||
|
|
||
| A - - - - - - B - | ||
| - * * * * * * * - | ||
| - * * * * * * * - | ||
| - * * * * * * * - | ||
| - * * * * * * * - | ||
| - * * * * * * * - | ||
| C * * * * * * D - | ||
| - - - - - - - - - | ||
|
|
||
| The sum of pixels in the area denoted by * is given by S = D + A - B - C. | ||
| """ | ||
| immutable IntegralArray{T, N, A} <: AbstractArray{T, N} | ||
| data::A | ||
| end | ||
|
|
||
| function IntegralArray(array::AbstractArray) | ||
| integral_array = Array{Images.accum(eltype(array))}(size(array)) | ||
| sd = coords_spatial(array) | ||
| cumsum!(integral_array, array, sd[1]) | ||
| for i = 2:length(sd) | ||
| cumsum!(integral_array, integral_array, sd[i]) | ||
| end | ||
| IntegralArray{eltype(array), ndims(array), typeof(array)}(integral_array) | ||
| end | ||
|
|
||
| linearindexing(A::IntegralArray) = Base.LinearFast() | ||
| size(A::IntegralArray) = size(A.data) | ||
| getindex(A::IntegralArray, i::Int...) = A.data[i...] | ||
|
|
||
| getindex(A::IntegralArray, ids::Tuple...) = getindex(A, ids[1]...) | ||
|
There was a problem hiding this comment. This is not commonly supported, what is the motivation? There was a problem hiding this comment. This was for the case where ClosedInterval is of N dimensions. If you do CartesianIndex(N) ± x, it returns a Tuple of ClosedIntervals. I could not find a way to splatter the return value (it always returns a Tuple). There was a problem hiding this comment. Oh, very interesting indeed. One could presumably force the user to write it as |
||
|
|
||
| function getindex(A::IntegralArray, i::ClosedInterval...) | ||
|
There was a problem hiding this comment. Best to declare this with two indexes, getindex{T,N}(A::IntegralArray{T,N}, i::Vararg{ClosedInterval,N})which will require that the number of indices matches the dimensionality of the array. (This is Julia PR #11242.) There was a problem hiding this comment. Probably also want to insist they are There was a problem hiding this comment. Yes, I forgot to add that. It cannot be a simple bilinear interpolation though as it is incorrect. I'll think of some other way to interpolate. Any ideas? There was a problem hiding this comment. Agreed this may not be entirely simple. Of course you should be able to calculate the actual pixel values from the integral image, and from that you should be able to calculate the actual sub-pixel integral image (where you assume the image to be piecewise constant, aka nearest-neighbor interpolation). I haven't sat down and tried to do this calculation myself yet, so I don't have a sense for how messy it will be. But the advantage is that it will automatically be consistent with the on-pixel integral image, which is a good thing. |
||
| _boxdiff(A, i[1].left, i[2].left, i[1].right, i[2].right) | ||
| end | ||
|
|
||
| function _boxdiff{T}(int_array::IntegralArray{T, 2}, tl_y::Integer, tl_x::Integer, br_y::Integer, br_x::Integer) | ||
| sum = int_array[br_y, br_x] | ||
| sum -= tl_x > 1 ? int_array[br_y, tl_x - 1] : zero(T) | ||
| sum -= tl_y > 1 ? int_array[tl_y - 1, br_x] : zero(T) | ||
| sum += tl_y > 1 && tl_x > 1 ? int_array[tl_y - 1, tl_x - 1] : zero(T) | ||
| sum | ||
| end | ||
|
|
||
| end # module | ||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1 @@ | ||
| IntervalSets |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,5 +1,52 @@ | ||
| module IntegralArraysTests | ||
|
|
||
| using IntegralArrays, Base.Test, IntervalSets | ||
|
|
||
| @testset "IntegralArrays" begin | ||
| a = zeros(10, 10) | ||
| int_array = IntegralArray(a) | ||
| @test all(int_array == a) | ||
|
|
||
| a = ones(10,10) | ||
| int_array = IntegralArray(a) | ||
| chk = Array(1:10) | ||
| @test all([vec(int_array[i, :]) == chk * i for i in 1:10]) | ||
|
|
||
| int_sum = int_array[1..5, 1..2] | ||
| @test int_sum == 10.0 | ||
| int_sum = int_array[3 ± 2, 1..2] | ||
| @test int_sum == 10.0 | ||
| int_sum = int_array[CartesianIndex((3, 3)) ± 2] | ||
| @test int_sum == 25.0 | ||
| int_sum = int_array[1..2, 1..5] | ||
| @test int_sum == 10.0 | ||
| int_sum = int_array[4..8, 4..8] | ||
| @test int_sum == 25.0 | ||
| int_sum = int_array[6 ± 2, 6 ± 2] | ||
| @test int_sum == 25.0 | ||
| int_sum = int_array[CartesianIndex((6, 6)) ± 2] | ||
| @test int_sum == 25.0 | ||
|
|
||
| a = Array(reshape(1:100, 10, 10)) | ||
| int_array = IntegralArray(a) | ||
| @test int_array[diagind(int_array)] == Array([1, 26, 108, 280, 575, 1026, 1666, 2528, 3645, 5050]) | ||
|
|
||
| int_sum = int_array[1..3, 1..3] | ||
| @test int_sum == 108 | ||
| int_sum = int_array[2 ± 1, 2 ± 1] | ||
| @test int_sum == 108 | ||
| int_sum = int_array[CartesianIndex((2, 2)) ± 1] | ||
| @test int_sum == 108 | ||
| int_sum = int_array[1..5, 1..2] | ||
| @test int_sum == 80 | ||
| int_sum = int_array[3 ± 2, 1..2] | ||
| @test int_sum == 80 | ||
| int_sum = int_array[4..8, 4..8] | ||
| @test int_sum == 1400 | ||
| int_sum = int_array[6 ± 2, 6 ± 2] | ||
| @test int_sum == 1400 | ||
| int_sum = int_array[CartesianIndex((6, 6)) ± 2] | ||
| @test int_sum == 1400 | ||
| end | ||
|
|
||
| end |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should we support indexing with integers at all? What use is there for the integral itself?
If we don't have to support Int, it might be good to have it throw a helpful error message (rather than just throw
MethodError).There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess my point is that folks might unthinkingly index with integers and assume it covers the range from the beginning to the indicated position.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Completely agree on this point. Will change it 😄