Skip to content

Commit ac4eb52

Browse files
emmtemmt
committed
Fix shifting and ambiguities in separable.jl 
1 parent fed2099 commit ac4eb52

File tree

2 files changed

+185
-191
lines changed

2 files changed

+185
-191
lines changed

‎src/separable.jl

+120-169
Original file line numberDiff line numberDiff line change
@@ -164,22 +164,15 @@ where the second dimension is omitted from the list of dimensions.
164164
The *local average* of the two-dimensional array `A` on a centered moving window of size
165165
11×11 can be computed as:
166166
167-
localfilter(A, :, +, (-5:5, -5:5))*(1/11)
167+
localfilter(A, :, +, (-5:5, -5:5)) ./ 11
168168
169169
See [`localfilter!`](@ref) for an in-place version of the method.
170170
171171
"""
172-
function localfilter(A::AbstractArray,
173-
dims::Dimensions,
174-
op::Function, args...)
175-
return localfilter(eltype(A), A, dims, op, args...)
176-
end
177-
178-
function localfilter(T::Type, A::AbstractArray,
179-
dims::Dimensions,
180-
op::Function, args...)
181-
return localfilter!(similar(A, T), A, dims, op, args...)
182-
end
172+
localfilter(A::AbstractArray, dims::Dimensions, op::Function, args...) =
173+
localfilter(eltype(A), A, dims, op, args...)
174+
localfilter(::Type{T}, A::AbstractArray, dims::Dimensions, op::Function, args...) where {T} =
175+
localfilter!(similar(A, T), A, dims, op, args...)
183176

184177
"""
185178
localfilter!([dst = A,] A, dims, op, [ord=FORWARD_FILTER,] rngs[, wrk])
@@ -200,27 +193,24 @@ structuring element of width 7 in every dimension can be obtained by:
200193
201194
""" localfilter!
202195

203-
# In-place operation.
204-
localfilter!(A::AbstractArray, dims::Dimensions, op::Function, args...) =
205-
localfilter!(A, A, dims, op, args...)
206-
207196
# Provide ordering.
208-
function localfilter!(dst::AbstractArray{T,N},
209-
A::AbstractArray{<:Any,N},
210-
dims::Dimensions,
211-
op::Function,
212-
rngs::Ranges,
213-
args...) where {T,N}
197+
function localfilter!(A::AbstractArray, dims::Dimensions, op::Function, rngs::Ranges, args...)
198+
return localfilter!(A, dims, op, FORWARD_FILTER, rngs, args...)
199+
end
200+
function localfilter!(dst::AbstractArray{<:Any,N}, A::AbstractArray{<:Any,N},
201+
dims::Dimensions, op::Function, rngs::Ranges, args...) where {N}
214202
return localfilter!(dst, A, dims, op, FORWARD_FILTER, rngs, args...)
215203
end
216204

205+
# In-place operation.
206+
function localfilter!(A::AbstractArray, dims::Dimensions, op::Function,
207+
ord::FilterOrdering, rngs::Ranges, args...)
208+
return localfilter!(A, A, dims, op, ord, rngs, args...)
209+
end
210+
217211
# Provide workspace.
218-
function localfilter!(dst::AbstractArray{T,N},
219-
A::AbstractArray{<:Any,N},
220-
dims::Dimensions,
221-
op::Function,
222-
ord::FilterOrdering,
223-
rngs::Ranges) where {T,N}
212+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, dims::Dimensions,
213+
op::Function, ord::FilterOrdering, rngs::Ranges) where {T,N}
224214
wrk = Array{T}(undef, workspace_length(A, dims, rngs))
225215
return localfilter!(dst, A, dims, op, ord, rngs, wrk)
226216
end
@@ -230,126 +220,102 @@ end
230220
# The filter is applied along all chosen dimensions. To reduce page memory faults, the
231221
# operation is performed out-of-place (unless dst and A are the same) for the first chosen
232222
# dimension and the operation is performed in-place for the other chosen dimensions.
233-
234-
# Versions for a single given filter range.
235-
236-
function localfilter!(dst::AbstractArray{T,N},
237-
A::AbstractArray{<:Any,N},
238-
dims::Dimensions,
239-
op::Function,
240-
ord::FilterOrdering,
241-
rng::Axis,
242-
wrk::Vector{T}) where {T,N}
243-
# small optimization: convert range once
244-
return localfilter!(dst, A, dims, op, FORWARD_FILTER,
245-
kernel_range(ord, rng), wrk)
246-
end
247-
248-
function localfilter!(dst::AbstractArray{T,N},
249-
A::AbstractArray{<:Any,N},
250-
::Colon,
251-
op::Function,
252-
ord::ForwardFilterOrdering,
253-
rng::AbstractUnitRange{Int},
223+
#
224+
# Apply filter along all dimensions (`dims` is a colon).
225+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, ::Colon,
226+
op::Function, ord::FilterOrdering, len::Integer,
254227
wrk::Vector{T}) where {T,N}
255-
localfilter!(dst, A, 1, op, ord, rng, wrk)
256-
for d in 2:N
257-
localfilter!(dst, d, op, ord, rng, wrk)
258-
end
259-
return dst
228+
return localfilter!(dst, A, :, op, ord, kernel_range(len), wrk)
260229
end
261-
262-
function localfilter!(dst::AbstractArray{T,N},
263-
A::AbstractArray{<:Any,N},
264-
dims::Union{AbstractVector{<:Integer},
265-
Tuple{Vararg{Integer}}},
266-
op::Function,
267-
ord::ForwardFilterOrdering,
268-
rng::AbstractUnitRange{Int},
230+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, ::Colon,
231+
op::Function, ord::FilterOrdering, rng::AbstractRange{<:Integer},
269232
wrk::Vector{T}) where {T,N}
270-
if length(dims) ≥ 1
271-
i = firstindex(dims)
272-
localfilter!(dst, A, dims[i], op, ord, rng, wrk)
273-
while i < lastindex(dims)
274-
i += 1
275-
localfilter!(dst, dims[i], op, ord, rng, wrk)
233+
if N ≥ 1
234+
rng = kernel_range(ord, rng) # optimization: convert once
235+
localfilter!(dst, A, 1, op, FORWARD_FILTER, rng, wrk)
236+
for dim in 2:N
237+
localfilter!(dst, dst, dim, op, FORWARD_FILTER, rng, wrk)
276238
end
277239
end
278240
return dst
279241
end
280-
281-
282-
# Versions for a multiple given filter ranges.
283-
284-
function localfilter!(dst::AbstractArray{T,N},
285-
A::AbstractArray{<:Any,N},
286-
::Colon,
287-
op::Function,
288-
ord::FilterOrdering,
289-
rngs::Union{AbstractVector{<:Axis},
290-
Tuple{Vararg{Axis}}},
242+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, ::Colon,
243+
op::Function, ord::FilterOrdering, rngs::Ranges,
291244
wrk::Vector{T}) where {T,N}
292245
length(rngs) == N || throw(DimensionMismatch(
293246
"there must be as many intervals as dimensions"))
294247
if N ≥ 1
295-
dim = 1
296-
off = firstindex(rngs) - dim
297-
localfilter!(dst, A, dim, op, ord, rngs[off+dim], wrk)
298-
while dim < N
299-
dim += 1
300-
localfilter!(dst, dim, op, ord, rngs[off+dim], wrk)
248+
localfilter!(dst, A, 1, op, ord, first(rngs), wrk)
249+
for (dim, rng) in enumerate(rngs)
250+
dim > 1 && localfilter!(dst, dst, dim, op, ord, rng, wrk)
301251
end
302252
end
303253
return dst
304254
end
305-
306-
function localfilter!(dst::AbstractArray{T,N},
307-
A::AbstractArray{<:Any,N},
308-
dims::Union{AbstractVector{<:Integer},
309-
Tuple{Vararg{Integer}}},
310-
op::Function,
311-
ord::FilterOrdering,
312-
rngs::Union{AbstractVector{<:Axis},
313-
Tuple{Vararg{Axis}}},
255+
#
256+
# Apply filter along a single given dimension (`dims` is an integer).
257+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, dim::Integer,
258+
op::Function, ord::FilterOrdering, len::Integer,
314259
wrk::Vector{T}) where {T,N}
315-
len = length(dims)
316-
length(rngs) == len || throw(DimensionMismatch(
317-
"list of dimensions and list of intervals must have the same length"))
318-
if len ≥ 1
319-
idx = firstindex(dims)
320-
off = firstindex(rngs) - idx
321-
localfilter!(dst, A, dims[idx], op, ord, rngs[off+idx], wrk)
322-
while idx < lastindex(dims)
323-
idx += 1
324-
localfilter!(dst, dims[idx], op, ord, rngs[off+idx], wrk)
260+
return localfilter!(dst, A, dim, op, ord, kernel_range(len), wrk)
261+
end
262+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, dim::Integer,
263+
op::Function, ord::FilterOrdering, rng::AbstractRange{<:Integer},
264+
wrk::Vector{T}) where {T,N}
265+
# NOTE This is the most basic version which is called by other versions.
266+
1 ≤ dim ≤ N || throw(ArgumentError("out of bounds dimension"))
267+
isempty(rng) && throw(ArgumentError("invalid filter size"))
268+
unsafe_localfilter!(dst, A, Val{Int(dim)}(), op, kernel_range(ord, rng), wrk)
269+
return dst
270+
end
271+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, dim::Integer,
272+
op::Function, ord::FilterOrdering, rngs::Ranges,
273+
wrk::Vector{T}) where {T,N}
274+
length(rngs) == 1 || throw(DimensionMismatch(
275+
"there must be as many intervals as dimensions"))
276+
return localfilter!(dst, A, dim, op, ord, first(rngs), wrk)
277+
end
278+
#
279+
# Apply filter along multiple given dimensions (`dims` is a list of integers).
280+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, dims::Dimensions,
281+
op::Function, ord::FilterOrdering, len::Integer,
282+
wrk::Vector{T}) where {T,N}
283+
return localfilter!(dst, A, dims, op, ord, kernel_range(len), wrk) # FIXME this one replaces all
284+
end
285+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, dims::Dimensions,
286+
op::Function, ord::FilterOrdering, rng::AbstractRange{<:Integer},
287+
wrk::Vector{T}) where {T,N}
288+
i, j = firstindex(dims), lastindex(dims)
289+
if i ≤ j
290+
rng = kernel_range(ord, rng) # optimization: convert once
291+
localfilter!(dst, A, @inbounds(dims[i]), op, FORWARD_FILTER, rng, wrk)
292+
for k in i+1:j
293+
localfilter!(dst, dst, @inbounds(dims[k]), op, FORWARD_FILTER, rng, wrk)
325294
end
326295
end
327296
return dst
328297
end
329-
330-
# Apply filter along a single given dimension. This is the most basic version which is
331-
# called by other versions.
332-
function localfilter!(dst::AbstractArray{T,N},
333-
A::AbstractArray{<:Any,N},
334-
dim::Integer, # dimension of interest
335-
op::Function,
336-
ord::FilterOrdering,
337-
rng::Axis,
298+
function localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N}, dims::Dimensions,
299+
op::Function, ord::FilterOrdering, rngs::Ranges,
338300
wrk::Vector{T}) where {T,N}
339-
1 ≤ dim ≤ N || throw(ArgumentError("out of bounds dimension"))
340-
isempty(rng) && throw(ArgumentError("invalid filter size"))
341-
unsafe_localfilter!(dst, A, Val(Int(dim)), op, kernel_range(ord, rng), wrk)
301+
length(dims) == length(rngs) || throw(DimensionMismatch(
302+
"list of dimensions and list of intervals must have the same length"))
303+
i, j = firstindex(dims), lastindex(dims)
304+
if i ≤ j
305+
off = firstindex(rngs) - i
306+
localfilter!(dst, A, @inbounds(dims[i]), op, ord, @inbounds(rngs[off+i]), wrk)
307+
for k in i+1:j
308+
localfilter!(dst, dst, @inbounds(dims[k]), op, ord, @inbounds(rngs[off+k]), wrk)
309+
end
310+
end
342311
return dst
343312
end
344313

345314
# This version is to break the type instability related to the variable dimension of
346315
# interest.
347-
function unsafe_localfilter!(dst::AbstractArray{T,N},
348-
A::AbstractArray{<:Any,N},
349-
::Val{D}, # dimension of interest
350-
op::Function,
351-
K::AbstractUnitRange{Int},
352-
wrk::Vector{T}) where {T,N,D}
316+
function unsafe_localfilter!(dst::AbstractArray{T,N}, A::AbstractArray{<:Any,N},
317+
::Val{D} #= dimension of interest =#, op::Function,
318+
K::AbstractUnitRange{Int}, wrk::Vector{T}) where {T,N,D}
353319
src_inds = axes(A)
354320
dst_inds = axes(dst)
355321
for d in 1:D-1
@@ -386,15 +352,15 @@ function unsafe_shiftarray!(dst::AbstractArray{<:Any,N},
386352
k::Int) where {N}
387353
@inbounds for i2 ∈ I2, i1 ∈ I1
388354
# To allow for in-place operation without using any temporaries, we walk along the
389-
# dimension in the forward direction if the shift is nonnegative and in the
390-
# reverse direction otherwise.
391-
if k ≥ 0
392-
@simd for i ∈ I
393-
dst[i1,i,i2] = src[i1,B(i+k),i2]
355+
# dimension in the reverse direction if the shift is positive and in the forward
356+
# direction otherwise.
357+
if k > 0
358+
@simd for i ∈ reverse(I)
359+
dst[i1,i,i2] = src[i1,B(i-k),i2]
394360
end
395361
else
396-
@simd for i ∈ reverse(I)
397-
dst[i1,i,i2] = src[i1,B(i+k),i2]
362+
@simd for i ∈ I
363+
dst[i1,i,i2] = src[i1,B(i-k),i2]
398364
end
399365
end
400366
end
@@ -580,41 +546,12 @@ no needs for a workspace array.
580546
workspace_length(n::Int, p::Int) = ifelse((n < 1)|(p ≤ 1), 0, n + 2*(p - 1))
581547
workspace_length(n::Integer, p::Integer) = workspace_length(Int(n), Int(p))
582548

583-
function workspace_length(A::AbstractArray, dims::Dimensions,
584-
rng::OrdinalRange{<:Integer,<:Integer})
585-
return workspace_length(A, dims, length(rng))
586-
end
587-
588-
function workspace_length(A::AbstractArray, dim::Integer, len::Integer)
589-
if 1 ≤ dim ≤ ndims(A)
590-
return workspace_length(length(axes(A, dim)), len)
591-
else
592-
return 0
593-
end
594-
end
595-
596-
function workspace_length(A::AbstractArray, ::Colon, len::Integer)
597-
if ndims(A) ≥ 1
598-
return workspace_length(reduce(max, map(length, axes(A))), len)
599-
else
600-
return 0
601-
end
602-
end
603-
604-
function workspace_length(A::AbstractArray,
605-
dims::Union{AbstractVector{<:Integer},
606-
Tuple{Vararg{Integer}}},
607-
len::Integer)
608-
result = 0
609-
for dim in dims
610-
result = max(result, workspace_length(A, dim, len))
611-
end
612-
return result
613-
end
614-
615-
function workspace_length(A::AbstractArray, ::Colon,
616-
rngs::Union{AbstractVector{<:Axis},
617-
Tuple{Vararg{Axis}}})
549+
# `dims` is a colon.
550+
workspace_length(A::AbstractArray, ::Colon, rng::AbstractRange{<:Integer}) =
551+
workspace_length(A, :, length(rng))
552+
workspace_length(A::AbstractArray, ::Colon, len::Integer) =
553+
ndims(A) ≥ 1 ? workspace_length(maximum(size(A)), len) : 0
554+
function workspace_length(A::AbstractArray, ::Colon, rngs::Ranges)
618555
result = 0
619556
if length(rngs) == ndims(A)
620557
for (dim, rng) in enumerate(rngs)
@@ -624,11 +561,25 @@ function workspace_length(A::AbstractArray, ::Colon,
624561
return result
625562
end
626563

627-
function workspace_length(A::AbstractArray,
628-
dims::Union{AbstractVector{<:Integer},
629-
Tuple{Vararg{Integer}}},
630-
rngs::Union{AbstractVector{<:Axis},
631-
Tuple{Vararg{Axis}}})
564+
# `dims` is a single dimension.
565+
workspace_length(A::AbstractArray, dim::Integer, rng::AbstractRange{<:Integer}) =
566+
workspace_length(A, dim, length(rng))
567+
workspace_length(A::AbstractArray, dim::Integer, len::Integer) =
568+
1 ≤ dim ≤ ndims(A) ? workspace_length(size(A, dim), len) : 0
569+
workspace_length(A::AbstractArray, dim::Integer, rngs::Ranges) =
570+
length(rngs) == 1 ? workspace_length(A, dim, first(rngs)) : 0
571+
572+
# `dims` is a list of dimensions.
573+
workspace_length(A::AbstractArray, dims::Dimensions, rng::AbstractRange{<:Integer}) =
574+
workspace_length(A, dims, length(rng))
575+
function workspace_length(A::AbstractArray, dims::Dimensions, len::Integer)
576+
result = 0
577+
for dim in dims
578+
result = max(result, workspace_length(A, dim, len))
579+
end
580+
return result
581+
end
582+
function workspace_length(A::AbstractArray, dims::Dimensions, rngs::Ranges)
632583
result = 0
633584
if length(rngs) == length(dims)
634585
for (dim, rng) in zip(dims, rngs)

0 commit comments

Comments
 (0)