Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix mapwindow for offset arrays #160

Merged
merged 6 commits into from
Jul 9, 2020
Merged

Fix mapwindow for offset arrays #160

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
77cb511
Fix mapwindow for offset arrays
yha Mar 23, 2020
ef313e0
Fix assuming merging of OffsetArrays#105 (IdOffsetRange preserve indi…
yha May 25, 2020
92a4aa8
Fix output axes for border=Inner(), test with anonymous functions
yha May 26, 2020
31c3312
update OffsetArrays dependency
johnnychen94 Jun 17, 2020
7663f48
Merge branch 'master' into mapwindow-offset
johnnychen94 Jun 17, 2020
5bab482
Clearer comment. Fix eltype in _indexof.
yha Jun 24, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "ImageFiltering"
uuid = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
author = ["Tim Holy <tim.holy@gmail.com>", "Jan Weidner <jw3126@gmail.com>"]
version = "0.6.13"
version = "0.6.14"

[deps]
CatIndices = "aafaddc9-749c-510e-ac4f-586e18779b91"
Expand Down Expand Up @@ -30,7 +30,7 @@ FFTW = "0.3, 1"
ImageCore = "0.8.1"
ImageMetadata = "0.9" # it doesn't actually depend on ImageMetadata, but if present we want arraydata to work
MappedArrays = "0.2"
OffsetArrays = "0.10, 0.11, 1"
OffsetArrays = "1.1"
Requires = "0.5, 1.0"
StaticArrays = "0.10, 0.11, 0.12"
TiledIteration = "0.2"
Expand Down
43 changes: 23 additions & 20 deletions src/mapwindow.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ using ..ImageFiltering: BorderSpecAny, Pad, Fill, Inner,
borderinstance, _interior, padindex, imfilter
using Base: Indices, tail
using Statistics
using OffsetArrays

export mapwindow, mapwindow!

Expand Down Expand Up @@ -72,14 +73,18 @@ function mapwindow(f, img, window; border="replicate",
resolve_imginds(indices))
end

# a unit range `r` having `r == axes(r)` which keeps its axes with
# broadcasting (e.g., `axes(r.+1) == axes(r)`), unlike Base.IdentityUnitRange
self_offset(r::AbstractUnitRange) = OffsetArrays.IdOffsetRange(1:length(r), first(r)-one(eltype(r)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not exactly sure I understand what you're trying to do. It doesn't quite do what it claims in the comment, for example

julia> r2 = OffsetArrays.IdOffsetRange(5:10, -3)
2:7

julia> r3 = self_offset(r2)
2:7

julia> first(r3)
2

julia> first(axes(r2, 1))
-2

What's the underlying goal? Would you like to create a range that preserves the values and indices of an original AbstractUnitRange but which behaves as described under broadcasting?

Or is this a situation where we can essentially ignore certain axes? For instance, we might imagine ignoring the axes of imginds, although alternatively they might specify the axes of the output.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suppose the comment is confusing because it uses the name r for the output of the function when it's also as the name of the function parameter. It should really say that s = self_offset(r) has the same values as r but has itself as indices: collect(s) == collect(r) and axes(s) == (s,) (in your example, axes(r3) == (2:7,) == (r3,)).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed the comment


function default_imginds(img, window, border)
axes(img)
end
function default_imginds(img, window, border::Inner)
imginds = axes(img)
win = resolve_window(window)
indind = _indices_of_interiour_indices(imginds, imginds, win)
map((II, r) -> r[II], indind, imginds)
indind = _indices_of_interiour_indices(imginds, imginds, win)::Tuple{Vararg{AbstractUnitRange}}
map(self_offset, indind)
end

function _mapwindow(f, img, window, border, imginds)
Expand Down Expand Up @@ -153,28 +158,26 @@ end
struct _IdentityTransformer <: _IndexTransformer end
@inline Base.getindex(t::_IdentityTransformer, inds) = inds

function _IndexTransformer(ranges)
stride = map(step, ranges)
offset1 = map(first, ranges)
offset = offset1 .- stride
function _IndexTransformer(from_ranges::NTuple{N,AbstractUnitRange}, to_ranges) where {N}
stride = map(step, to_ranges)
offset = first.(to_ranges) .- first.(from_ranges) .* stride
_AffineTransformer(offset, stride)
end

function _IndexTransformer(ranges::NTuple{N,Base.OneTo}) where {N}
_IdentityTransformer()
function _IndexTransformer(from_ranges::NTuple{N,AbstractUnitRange}, to_ranges::NTuple{N,AbstractUnitRange}) where {N}
offset = first.(to_ranges) .- first.(from_ranges)
_OffsetTransformer(offset)
end

function _IndexTransformer(ranges::NTuple{N,AbstractUnitRange}) where {N}
offset1 = map(first, ranges)
offset = offset1 .- 1
_OffsetTransformer(offset)
function _IndexTransformer(::NTuple{N,Base.OneTo}, ::NTuple{N,Base.OneTo}) where {N}
_IdentityTransformer()
end

compute_output_range(r::AbstractUnitRange) = r
compute_output_range(r::AbstractRange) = Base.OneTo(length(r))

function compute_output_indices(imginds)
ranges = map(compute_output_range, imginds)
ranges = map(i->compute_output_range(axes(i,1)), imginds)
# Base.similar does not like if some but not all ranges are Base.OneTo
homogenize(ranges)
end
Expand All @@ -190,14 +193,14 @@ function _intersectionindices(full::AbstractUnitRange, r::AbstractRange)
else
ret = _indexof(r,first(r_sub)):_indexof(r,last(r_sub))
end
@assert intersect(full, r) == r[ret]
@assert r_sub == r[ret] || isempty(r_sub) && isempty(ret)
ret
end

function _indexof(r::AbstractRange, x)
T = eltype(r)
yha marked this conversation as resolved.
Show resolved Hide resolved
@assert x ∈ r
i = one(T) + T((x - first(r)) / step(r))
i = T(firstindex(r) + (x - first(r)) / step(r))
@assert r[i] == x
i
end
Expand Down Expand Up @@ -258,13 +261,13 @@ function mapwindow_kernel!(f,

@assert map(length, imginds) == map(length, axes(out))

indind_full = map(r -> Base.OneTo(length(r)), imginds)
indind_full = map(r -> axes(r,1), imginds)
indind_inner = _indices_of_interiour_indices(axes(img), imginds, window)
Rindind_full = CartesianIndices(indind_full)
Rindind_inner = CartesianIndices(indind_inner)

outindtrafo = _IndexTransformer(axes(out))
imgindtrafo = _IndexTransformer(imginds)
outindtrafo = _IndexTransformer(indind_full, axes(out))
imgindtrafo = _IndexTransformer(indind_full, imginds)

buf, bufrs = allocate_buffer(f, img, window)
Rbuf = CartesianIndices(size(buf))
Expand Down Expand Up @@ -375,7 +378,7 @@ extrema_filter(A::AbstractArray, window) = error("`window` must have the same nu

extrema_filter(A::AbstractArray{T,N}, window::Integer) where {T,N} = extrema_filter(A, ntuple(d->window, Val{N}))

function _extrema_filter!(A::Array, w1, w...)
function _extrema_filter!(A::AbstractArray, w1, w...)
if w1 > 1
a = first(A)
if w1 <= 20
Expand All @@ -396,7 +399,7 @@ function _extrema_filter!(A::Array, w1, w...)
return A
end
end
_extrema_filter!(A::Array) = A
_extrema_filter!(A::AbstractArray) = A

# Extrema-filtering along "columns" (dimension 1). This implements Lemire
# Algorithm 1, with the following modifications:
Expand Down
31 changes: 30 additions & 1 deletion test/mapwindow.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,33 @@ using ImageFiltering: IdentityUnitRange
Amin = groundtruth(min, A, w)
@test minval == Amin
end
# offsets
@testset "offsets" begin
n = 5
arrays = [rand(n), rand(n,n), rand(n,n,n)]
@testset "offsets ($f, $offset, $window, $dim)" for f in (extrema,maximum,minimum),
offset in -5:5,
window in [1:2:9; [0:2,-2:0]],
# broken: "reflect", NA(), NoPad()
border in ["replicate", "symmetric", Fill(randn()), Inner()],
(dim,a) in enumerate(arrays)
offsets = ntuple(_->offset,dim)
windows = ntuple(_->window,dim)
winlen = window isa Number ? window : length(window)
wrapped_f = x->f(x)
for g in (f, wrapped_f)
mw(a) = mapwindow(g, a, windows, border=border)

if border == Inner() && winlen > n
@test_throws DimensionMismatch mw(a)
@test_throws DimensionMismatch mw(OffsetArray(a,offsets))
else
@test OffsetArray(mw(a),offsets) == mw(OffsetArray(a,offsets))
end
end
end
end
johnnychen94 marked this conversation as resolved.
Show resolved Hide resolved


# median
for f in (median, median!)
Expand Down Expand Up @@ -115,7 +142,7 @@ using ImageFiltering: IdentityUnitRange
@test expected == @inferred mapwindow!(f,out,img,window,indices=imginds)
end
for (inds, kw) ∈ [((Base.OneTo(3),), (border="replicate", indices=2:2:7)),
((2:7,), (indices=2:7,)),
(axes(2:7), (indices=2:7,)),
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do not have an opinion on whether this is a bug fix or breaking change.

((Base.OneTo(10),), ())
]
@test inds == axes(mapwindow(mean, randn(10), (3,); kw...))
Expand All @@ -130,6 +157,8 @@ using ImageFiltering: IdentityUnitRange
@test mapwindow(first, img_48, (0:2,),
border=Inner(),
indices=inds_48) == img_48[inds_48]
res_shifted = mapwindow(minimum, img_48, -2:0, border=Inner())
@test res_shifted == OffsetArray(img_48[1:8], 3:10)

@testset "desugaring window argument #58" begin
img58 = rand(10)
Expand Down