From 6e8901b3b03210905081bd9826562c299dfd068e Mon Sep 17 00:00:00 2001 From: platawiec Date: Wed, 4 Jul 2018 02:30:21 -0400 Subject: [PATCH] add n-dimensional unwrap (#198) --- src/DSP.jl | 3 +- src/Filters/Filters.jl | 1 + src/Filters/response.jl | 2 +- src/unwrap.jl | 378 ++++++++++++++++++++++++++++++++++++++++ src/util.jl | 65 +------ test/runtests.jl | 2 +- test/unwrap.jl | 140 +++++++++++++++ test/util.jl | 40 ----- 8 files changed, 524 insertions(+), 107 deletions(-) create mode 100644 src/unwrap.jl create mode 100644 test/unwrap.jl diff --git a/src/DSP.jl b/src/DSP.jl index 068d0d74d..9c42e5e4a 100644 --- a/src/DSP.jl +++ b/src/DSP.jl @@ -46,6 +46,7 @@ if VERSION >= v"0.7.0-DEV.602" end include("util.jl") +include("unwrap.jl") include("windows.jl") include("periodograms.jl") include("Filters/Filters.jl") @@ -53,7 +54,7 @@ include("lpc.jl") include("estimation.jl") using Reexport -@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Estimation +@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Unwrap, .Estimation include("deprecated.jl") end diff --git a/src/Filters/Filters.jl b/src/Filters/Filters.jl index 9c4cbad1b..5f3731ecc 100644 --- a/src/Filters/Filters.jl +++ b/src/Filters/Filters.jl @@ -1,5 +1,6 @@ module Filters using ..DSP: @importffts, mul!, rmul! +using ..Unwrap using Polynomials, ..Util import Base: * using Compat: copyto!, undef diff --git a/src/Filters/response.jl b/src/Filters/response.jl index dd9138bec..08189aff6 100644 --- a/src/Filters/response.jl +++ b/src/Filters/response.jl @@ -54,7 +54,7 @@ or frequencies `w` in radians/sample. """ function phasez(filter::FilterCoefficients, w = Compat.range(0, stop=π, length=250)) h = freqz(filter, w) - unwrap(-atan2.(imag(h), real(h))) + unwrap(-atan2.(imag(h), real(h)); dims=ndims(h)) end diff --git a/src/unwrap.jl b/src/unwrap.jl new file mode 100644 index 000000000..2cf1acfb4 --- /dev/null +++ b/src/unwrap.jl @@ -0,0 +1,378 @@ +module Unwrap + +using Compat.Random + +import Compat +import Compat: selectdim +import Compat: CartesianIndices +import Compat: Nothing +import Compat.Random: GLOBAL_RNG + +export unwrap, unwrap! + +""" + unwrap!(m; kwargs...) + +In-place version of [`unwrap`](@ref). +""" +function unwrap!(m::AbstractArray{T,N}; dims=nothing, kwargs...) where {T, N} + if dims === nothing && N != 1 + Base.depwarn("`unwrap!(m::AbstractArray)` is deprecated, use `unwrap!(m, dims=ndims(m))` instead", :unwrap!) + dims = N + end + unwrap!(m, m; dims=dims, kwargs...) +end + +""" + unwrap!(y, m; kwargs...) + +Unwrap `m` storing the result in `y`, see [`unwrap`](@ref). +""" +function unwrap!(y::AbstractArray{T,N}, m::AbstractArray{T,N}; dims=nothing, range=2T(pi), kwargs...) where {T, N} + if dims === nothing + if N != 1 + throw(ArgumentError("`unwrap!`: required keyword parameter dims missing")) + end + dims = 1 + end + if dims isa Integer + Compat.accumulate!(unwrap_kernel(range), y, m, dims=dims) + elseif dims == 1:N + unwrap_nd!(y, m; range=range, kwargs...) + else + throw(ArgumentError("`unwrap!`: Invalid dims specified: $dims")) + end + return y +end + +unwrap_kernel(range) = (x, y) -> y - round((y - x) / range) * range + +""" + unwrap(m; kwargs...) + + +Assumes `m` to be a sequence of values that has been wrapped to be inside the +given `range` (centered around zero), and undoes the wrapping by identifying +discontinuities. If a single dimension is passed to `dims`, then `m` is assumed +to have wrapping discontinuities only along that dimension. If a range of +dimensions, as in `1:ndims(m)`, is passed to `dims`, then `m` is assumed to have +wrapping discontinuities across all `ndims(m)` dimensions. + +A common usage for unwrapping across a singleton dimension is for a phase +measurement over time, such as when +comparing successive frames of a short-time-fourier-transform, as +each frame is wrapped to stay within (-pi, pi]. + +A common usage for unwrapping across multiple dimensions is for a phase +measurement of a scene, such as when retrieving the phase information of +of an image, as each pixel is wrapped to stay within (-pi, pi]. + +# Arguments +- `m::AbstractArray{T, N}`: Array to unwrap. +- `dims=nothing`: Dimensions along which to unwrap. If `dims` is an integer, then + `unwrap` is called on that dimension. If `dims=1:ndims(m)`, then `m` is unwrapped + across all dimensions. +- `range=2pi`: Range of wrapped array. +- `circular_dims=(false, ...)`: When an element of this tuple is `true`, the + unwrapping process will consider the edges along the corresponding axis + of the array to be connected. +- `rng=GLOBAL_RNG`: Unwrapping of arrays with dimension > 1 uses a random + initialization. A user can pass their own RNG through this argument. +""" +function unwrap(m::AbstractArray{T,N}; dims=nothing, kwargs...) where {T, N} + if dims === nothing && N != 1 + Base.depwarn("`unwrap(m::AbstractArray)` is deprecated, use `unwrap(m, dims=ndims(m))` instead", :unwrap) + dims = ndims(m) + end + unwrap!(similar(m), m; dims=dims, kwargs...) +end + +@deprecate(unwrap(m::AbstractArray, dim::Integer; range::Number=2pi), + unwrap(m, dims=dim, range=range)) + +#= Algorithm based off of + M. A. Herráez, D. R. Burton, M. J. Lalor, and M. A. Gdeisat, + "Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path" + `Applied Optics, Vol. 41, Issue 35, pp. 7437-7444 (2002) ` + and + H. Abdul-Rahman, M. Gdeisat, D. Burton, M. Lalor, + "Fast three-dimensional phase-unwrapping algorithm based on sorting by reliability following a non-continuous path", + `Proc. SPIE 5856, Optical Measurement Systems for Industrial Inspection IV, 32 (2005) ` + Code inspired by Scipy's implementation, which is under BSD license. +=# + +mutable struct Pixel{T} + periods::Int + val::T + reliability::Float64 + groupsize::Int + head::Pixel{T} + last::Pixel{T} + next::Union{Nothing, Pixel{T}} + function Pixel{T}(periods, val, rel, gs) where T + pixel = new(periods, val, rel, gs) + pixel.head = pixel + pixel.last = pixel + pixel.next = nothing + return pixel + end +end +Pixel(v, rng) = Pixel{typeof(v)}(0, v, rand(rng), 1) +@inline Base.length(p::Pixel) = p.head.groupsize + +struct Edge{N} + reliability::Float64 + periods::Int + pixel_1::CartesianIndex{N} + pixel_2::CartesianIndex{N} +end +function Edge{N}(pixel_image::AbstractArray, ind1::CartesianIndex{N}, ind2::CartesianIndex{N}, range) where N + @inbounds rel = pixel_image[ind1].reliability + pixel_image[ind2].reliability + @inbounds periods = find_period(pixel_image[ind1].val, pixel_image[ind2].val, range) + return Edge{N}(rel, periods, ind1, ind2) +end +@inline Base.isless(e1::Edge, e2::Edge) = isless(e1.reliability, e2.reliability) + +function unwrap_nd!(dest::AbstractArray{T, N}, + src::AbstractArray{T, N}; + range::Number=2*convert(T, pi), + circular_dims::NTuple{N, Bool}=tuple(fill(false, N)...), + rng::AbstractRNG=GLOBAL_RNG) where {T, N} + + range_T = convert(T, range) + + pixel_image = init_pixels(src, rng) + calculate_reliability(pixel_image, circular_dims, range_T) + edges = Edge{N}[] + num_edges = _predict_num_edges(size(src), circular_dims) + sizehint!(edges, num_edges) + for idx_dim=1:N + populate_edges!(edges, pixel_image, idx_dim, circular_dims[idx_dim], range_T) + end + + sort!(edges, alg=MergeSort) + gather_pixels!(pixel_image, edges) + unwrap_image!(dest, pixel_image, range_T) + + return dest +end + +function _predict_num_edges(size_img, circular_dims) + num_edges = 0 + for (size_dim, wrap_dim) in zip(size_img, circular_dims) + num_edges += prod(size_img) * (size_dim-1) ÷ size_dim + wrap_dim * prod(size_img) ÷ size_dim + end + return num_edges +end + +# function to broadcast +function init_pixels(wrapped_image::AbstractArray{T, N}, rng) where {T, N} + pixel_image = similar(wrapped_image, Pixel{T}) + @Threads.threads for i in eachindex(wrapped_image) + @inbounds pixel_image[i] = Pixel(wrapped_image[i], rng) + end + return pixel_image +end + +function gather_pixels!(pixel_image, edges) + for edge in edges + @inbounds p1 = pixel_image[edge.pixel_1] + @inbounds p2 = pixel_image[edge.pixel_2] + merge_groups!(edge, p1, p2) + end +end + +function unwrap_image!(dest, pixel_image, range) + @Threads.threads for i in eachindex(dest) + @inbounds dest[i] = range * pixel_image[i].periods + pixel_image[i].val + end +end + +function wrap_val(val, range) + wrapped_val = val + wrapped_val += ifelse(val > range/2, -range, zero(val)) + wrapped_val += ifelse(val < -range/2, range, zero(val)) + return wrapped_val +end + +function find_period(val_left, val_right, range) + difference = val_left - val_right + period = 0 + period += ifelse(difference > range/2, -1, 0) + period += ifelse(difference < -range/2, 1, 0) + return period +end + +function merge_groups!(edge, pixel_1, pixel_2) + if is_differentgroup(pixel_1, pixel_2) + # pixel 2 is alone in group + if is_pixelalone(pixel_2) + merge_pixels!(pixel_1, pixel_2, -edge.periods) + elseif is_pixelalone(pixel_1) + merge_pixels!(pixel_2, pixel_1, edge.periods) + else + if is_bigger(pixel_1, pixel_2) + merge_into_group!(pixel_1, pixel_2, -edge.periods) + else + merge_into_group!(pixel_2, pixel_1, edge.periods) + end + end + end +end + +@inline function is_differentgroup(p1::Pixel, p2::Pixel) + return p1.head !== p2.head +end +@inline function is_pixelalone(pixel::Pixel) + return pixel.head === pixel.last +end +@inline function is_bigger(p1::Pixel, p2::Pixel) + return length(p1) ≥ length(p2) +end + +function merge_pixels!(pixel_base::Pixel, pixel_target::Pixel, periods) + pixel_base.head.groupsize += pixel_target.head.groupsize + pixel_base.head.last.next = pixel_target.head + pixel_base.head.last = pixel_target.head.last + pixel_target.head = pixel_base.head + pixel_target.periods = pixel_base.periods + periods +end + +function merge_into_group!(pixel_base::Pixel, pixel_target::Pixel, periods) + add_periods = pixel_base.periods + periods - pixel_target.periods + pixel = pixel_target.head + while pixel ≠ nothing + # merge all pixels in pixel_target's group to pixel_base's group + if pixel !== pixel_target + pixel.periods += add_periods + pixel.head = pixel_base.head + end + pixel = pixel.next + end + # assign pixel_target to pixel_base's group last + merge_pixels!(pixel_base, pixel_target, periods) +end + +function populate_edges!(edges, pixel_image::Array{T, N}, dim, connected, range) where {T, N} + size_img = collect(size(pixel_image)) + size_img[dim] -= 1 + idx_step = fill(0, N) + idx_step[dim] += 1 + idx_step_cart = CartesianIndex{N}(idx_step...) + idx_size = CartesianIndex{N}(size_img...) + for i in CartesianIndices(idx_size) + push!(edges, Edge{N}(pixel_image, i, i+idx_step_cart, range)) + end + if connected + idx_step = fill(0, N) + idx_step[dim] = -size_img[dim] + idx_step_cart = CartesianIndex{N}(idx_step...) + edge_begin = fill(1, N) + edge_begin[dim] = size(pixel_image)[dim] + edge_begin_cart = CartesianIndex{N}(edge_begin...) + for i in CartesianIndices(ntuple(dim_idx -> edge_begin_cart[dim_idx]:size(pixel_image, dim_idx), N)) + push!(edges, Edge{N}(pixel_image, i, i+idx_step_cart, range)) + end + end +end + +function calculate_reliability(pixel_image::AbstractArray{T, N}, circular_dims, range) where {T, N} + # get the shifted pixel indices in CartesinanIndex form + # This gets all the nearest neighbors (CartesionIndex{N}() = one(CartesianIndex{N})) + pixel_shifts = CartesianIndices(ntuple(i -> -1:1, N)) + size_img = size(pixel_image) + # inner loop + for i in CartesianIndices(ntuple(dim -> 2:(size(pixel_image, dim)-1), N)) + @inbounds pixel_image[i].reliability = calculate_pixel_reliability(pixel_image, i, pixel_shifts, range) + end + + for (idx_dim, connected) in enumerate(circular_dims) + if connected + # first border + pixel_shifts_border = copy(pixel_shifts) + for (idx_ps, ps) in enumerate(pixel_shifts_border) + # if the pixel shift goes out of bounds, we make the shift wrap + if ps[idx_dim] == 1 + new_ps = fill(0, N) + new_ps[idx_dim] = -size_img[idx_dim]+1 + pixel_shifts_border[idx_ps] = CartesianIndex{N}(new_ps...) + end + end + border_range = get_border_range(size_img, idx_dim, size_img[idx_dim]) + for i in CartesianIndices(border_range) + @inbounds pixel_image[i].reliability = calculate_pixel_reliability(pixel_image, i, pixel_shifts_border, range) + end + # second border + pixel_shifts_border = Compat.copyto!(pixel_shifts_border, pixel_shifts) + for (idx_ps, ps) in enumerate(pixel_shifts_border) + # if the pixel shift goes out of bounds, we make the shift wrap, this time to the other side + if ps[idx_dim] == -1 + new_ps = fill(0, N) + new_ps[idx_dim] = size_img[idx_dim]-1 + pixel_shifts_border[idx_ps] = CartesianIndex{N}(new_ps...) + end + end + border_range = get_border_range(size_img, idx_dim, 1) + for i in CartesianIndices(border_range) + @inbounds pixel_image[i].reliability = calculate_pixel_reliability(pixel_image, i, pixel_shifts_border, range) + end + end + end +end + +function get_border_range(size_img::NTuple{N, T}, border_dim, border_idx) where {N, T} + border_range = [2:(size_img[dim]-1) for dim=1:N] + border_range[border_dim] = border_idx:border_idx + return tuple(border_range...) +end + +function calculate_pixel_reliability(pixel_image::AbstractArray{Pixel{T}, N}, pixel_index, pixel_shifts, range) where {T, N} + sum_val = zero(T) + for pixel_shift in pixel_shifts + @inbounds sum_val += wrap_val(pixel_image[pixel_index+pixel_shift].val - pixel_image[pixel_index].val, range)^2 + end + return sum_val +end + +# specialized pixel reliability calculations for different N +function calculate_pixel_reliability(pixel_image::AbstractArray{Pixel{T}, 2}, pixel_index, pixel_shifts, range) where T + @inbounds D1 = wrap_val(pixel_image[pixel_index+pixel_shifts[2]].val - pixel_image[pixel_index].val, range) + @inbounds D2 = wrap_val(pixel_image[pixel_index+pixel_shifts[4]].val - pixel_image[pixel_index].val, range) + @inbounds H = wrap_val(pixel_image[pixel_index+pixel_shifts[6]].val - pixel_image[pixel_index].val, range) + @inbounds V = wrap_val(pixel_image[pixel_index+pixel_shifts[8]].val - pixel_image[pixel_index].val, range) + return H*H + V*V + D1*D1 + D2*D2 +end + +function calculate_pixel_reliability(pixel_image::AbstractArray{Pixel{T}, 3}, pixel_index, pixel_shifts, range) where T + sum_val = zero(T) + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[1]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[2]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[3]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[4]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[5]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[6]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[7]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[8]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[9]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[10]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[11]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[12]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[13]].val - pixel_image[pixel_index].val, range))^2 + # pixel_shifts[14] is null shift + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[15]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[16]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[17]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[18]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[19]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[20]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[21]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[22]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[23]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[24]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[25]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[26]].val - pixel_image[pixel_index].val, range))^2 + @inbounds sum_val += (wrap_val(pixel_image[pixel_index+pixel_shifts[27]].val - pixel_image[pixel_index].val, range))^2 + return sum_val +end + +end diff --git a/src/util.jl b/src/util.jl index 232088268..c5a9d0208 100644 --- a/src/util.jl +++ b/src/util.jl @@ -6,9 +6,7 @@ using Compat: copyto!, ComplexF64, selectdim, undef import Compat.LinearAlgebra.BLAS @importffts -export unwrap!, - unwrap, - hilbert, +export hilbert, Frequencies, fftintype, fftouttype, @@ -29,67 +27,6 @@ export unwrap!, polyfit, shiftin! -""" - unwrap!(m; dims=nothing, range=2pi) - -In-place version of [`unwrap`](@ref). -""" -function unwrap!(m::AbstractArray{<:Any,N}; dims=nothing, range=2pi) where N - if dims === nothing && N != 1 - Base.depwarn("`unwrap!(m::AbstractArray)` is deprecated, use `unwrap!(m, dims=ndims(m))` instead", :unwrap!) - dims = N - end - unwrap!(m, m; dims=dims, range=range) -end - -""" - unwrap!(y, m; dims=nothing, range=2pi) - -Unwrap `m` storing the result in `y`, see [`unwrap`](@ref). -""" -function unwrap!(y::AbstractArray{<:Any,N}, m::AbstractArray{<:Any,N}; dims=nothing, range::Number=2pi) where {N} - if dims === nothing - if N != 1 - throw(ArgumentError("`unwrap!`: required keyword parameter dims missing")) - end - dims = 1 - end - if dims isa Integer - Compat.accumulate!(unwrap_kernel(range), y, m, dims=dims) - else - # this is where #198 can hook in - throw(ArgumentError("`unwrap!`: Invalid dims specified: $dims")) - end - return y -end - -unwrap_kernel(range=2π) = (x, y) -> y - round((y - x) / range) * range - -@deprecate(unwrap!(m::AbstractArray, dim::Integer; range::Number=2pi), - unwrap!(m, dims=dim, range=range)) - - -""" - unwrap(m; dims=nothing, range=2pi) - -Assumes m (along the single dimension dims) to be sequences of values that -have been wrapped to be inside the given range (centered around -zero), and undoes the wrapping by identifying discontinuities. - -A common usage is for a phase measurement over time, such as when -comparing successive frames of a short-time-fourier-transform, as -each frame is wrapped to stay within (-pi, pi]. -""" -function unwrap(m::AbstractArray{<:Any,N}; dims=nothing, range=2pi) where {N} - if dims === nothing && N != 1 - Base.depwarn("`unwrap(m::AbstractArray)` is deprecated, use `unwrap(m, dims=ndims(m))` instead", :unwrap) - dims = ndims(m) - end - unwrap!(similar(m), m; dims=dims, range=range) -end - -@deprecate(unwrap(m::AbstractArray, dim::Integer; range::Number=2pi), - unwrap(m, dims=dim, range=range)) function hilbert(x::StridedVector{T}) where T<:FFTW.fftwReal # Return the Hilbert transform of x (a real signal). diff --git a/test/runtests.jl b/test/runtests.jl index 92624dcf4..2d88e2f5d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,7 @@ using Compat, DSP, AbstractFFTs, FFTW, Compat.Test testfiles = [ "dsp.jl", "util.jl", "windows.jl", "filter_conversion.jl", "filter_design.jl", "filter_response.jl", "filt.jl", "filt_stream.jl", - "periodograms.jl", "resample.jl", "lpc.jl", "estimation.jl"] + "periodograms.jl", "resample.jl", "lpc.jl", "estimation.jl", "unwrap.jl"] for testfile in testfiles eval(@__MODULE__, :(@testset $testfile begin include($testfile) end)) diff --git a/test/unwrap.jl b/test/unwrap.jl new file mode 100644 index 000000000..c99f34584 --- /dev/null +++ b/test/unwrap.jl @@ -0,0 +1,140 @@ +using DSP, Compat, Compat.Test + +@testset "Unwrap 1D" begin + @test unwrap([0.1, 0.2, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] + @test unwrap([0.1, 0.2 + 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] + @test unwrap([0.1, 0.2 - 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] + @test unwrap([0.1, 0.2 - 2pi, 0.3 - 2pi, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] + @test unwrap([0.1 + 2pi, 0.2, 0.3, 0.4]) ≈ [0.1 + 2pi, 0.2 + 2pi, 0.3 + 2pi, 0.4 + 2pi] + @test unwrap([0.1, 0.2 + 6pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] + + test_v = [0.1, 0.2, 0.3 + 2pi, 0.4] + res_v = unwrap(test_v) + @test test_v ≈ [0.1, 0.2, 0.3 + 2pi, 0.4] + res_v .= 0 + unwrap!(res_v, test_v) + @test res_v ≈ [0.1, 0.2, 0.3, 0.4] + @test test_v ≈ [0.1, 0.2, 0.3 + 2pi, 0.4] + unwrap!(test_v) + @test test_v ≈ [0.1, 0.2, 0.3, 0.4] + + # test unwrapping within multi-dimensional array + wrapped = [0.1, 0.2 + 2pi, 0.3, 0.4] + unwrapped = [0.1, 0.2, 0.3, 0.4] + wrapped = hcat(wrapped, wrapped) + unwrapped = hcat(unwrapped, unwrapped) + @test unwrap(wrapped, dims=2) ≈ wrapped + @test unwrap(wrapped, dims=1) ≈ unwrapped + @test unwrap!(copy(wrapped), dims=2) ≈ wrapped + @test unwrap!(copy(wrapped), dims=1) ≈ unwrapped + + # this should eventually default to the multi-dimensional case + @test_throws ArgumentError unwrap!(similar(wrapped), wrapped) + + # test unwrapping with other ranges + unwrapped = [1.0:100;] + wrapped = Float64[i % 10 for i in unwrapped] + @test unwrap(wrapped, range=10) ≈ unwrapped + + # test generically typed unwrapping + types = (Float32, Float64, BigFloat) + for T in types + srand(1234) + A_unwrapped = collect(Compat.range(0, stop=4convert(T, π), length=10)) + A_wrapped = A_unwrapped .% (2convert(T, π)) + + @test unwrap(A_wrapped) ≈ A_unwrapped + unwrap!(A_wrapped) + @test A_wrapped ≈ A_unwrapped + + A_unwrapped_range = collect(Compat.range(0, stop=4, length=10)) + test_range = convert(T, 2) + A_wrapped_range = A_unwrapped_range .% test_range + @test unwrap(A_wrapped_range; range=test_range) ≈ A_unwrapped_range + end +end + +# tests for multi-dimensional unwrapping +@testset "Unwrap 2D" begin + types = (Float32, Float64, BigFloat) + for T in types + srand(1234) + v_unwrapped = collect(Compat.range(0, stop=4convert(T, π), length=7)) + A_unwrapped = v_unwrapped .+ v_unwrapped' + A_wrapped = A_unwrapped .% (2convert(T, π)) + + test_unwrapped = unwrap(A_wrapped, dims=1:2) + d = first(A_unwrapped) - first(test_unwrapped) + @test (test_unwrapped .+ d) ≈ A_unwrapped + unwrap!(A_wrapped, dims=1:2) + d = first(A_unwrapped) - first(A_wrapped) + @test (A_wrapped .+ d) ≈ A_unwrapped + + v_unwrapped_range = collect(Compat.range(0, stop=4, length=7)) + A_unwrapped_range = v_unwrapped_range .+ v_unwrapped_range' + test_range = convert(T, 2) + A_wrapped_range = A_unwrapped_range .% test_range + + test_unwrapped_range = unwrap(A_wrapped_range, dims=1:2; range=test_range) + d = first(A_unwrapped_range) - first(test_unwrapped_range) + @test (test_unwrapped_range .+ d) ≈ A_unwrapped_range + + # Test circular_dims + # after unwrapping, pixels at borders should be equal to corresponding pixels + # on other side + circular_dims = (true, true) + wa_vec = Compat.range(0, stop=4convert(T, π), length=10) + wa_uw = wa_vec .+ zeros(10)' + # make periodic + wa_uw[end, :] = wa_uw[1, :] + wa_w = wa_uw .% (2π) + wa_test = unwrap(wa_w, dims=1:2, circular_dims=circular_dims, rng=MersenneTwister(0)) + # with wrap-around, the borders should be equal, but for this problem the + # image may not be recovered exactly + @test wa_test[:, 1] ≈ wa_test[:, end] + @test wa_test[end, :] ≈ wa_test[1, :] + # In this case, calling unwrap w/o circular_dims does not recover the borders + wa_test_nowa = unwrap(wa_w, dims=1:2) + @test !(wa_test_nowa[end, :] ≈ wa_test_nowa[1, :]) + + end +end + +@testset "Unwrap 3D" begin + types = (Float32, Float64, BigFloat) + f(x, y, z) = 0.1x^2 - 2y + 2z + f_wraparound2(x, y, z) = 5*sin(x) + 2*cos(y) + z + f_wraparound3(x, y, z) = 5*sin(x) + 2*cos(y) - 4*cos(z) + for T in types + grid = Compat.range(zero(T), stop=2convert(T, π), length=11) + f_uw = f.(grid, grid', reshape(grid, 1, 1, :)) + f_wr = f_uw .% (2convert(T, π)) + uw_test = unwrap(f_wr, dims=1:3) + offset = first(f_uw) - first(uw_test) + @test (uw_test.+offset) ≈ f_uw rtol=eps(T) #oop, nowrap + # test in-place version + unwrap!(f_wr, dims=1:3) + offset = first(f_uw) - first(f_wr) + @test (f_wr.+offset) ≈ f_uw rtol=eps(T) #ip, nowrap + + f_uw = f_wraparound2.(grid, grid', reshape(grid, 1, 1, :)) + f_wr = f_uw .% (2convert(T, π)) + uw_test = unwrap(f_wr, dims=1:3) + offset = first(f_uw) - first(uw_test) + @test (uw_test.+offset) ≈ f_uw #oop, 2wrap + # test in-place version + unwrap!(f_wr, dims=1:3, circular_dims=(true, true, false)) + offset = first(f_uw) - first(f_wr) + @test (f_wr.+offset) ≈ f_uw #ip, 2wrap + + f_uw = f_wraparound3.(grid, grid', reshape(grid, 1, 1, :)) + f_wr = f_uw .% (2convert(T, π)) + uw_test = unwrap(f_wr, dims=1:3, circular_dims=(true, true, true)) + offset = first(f_uw) - first(uw_test) + @test (uw_test.+offset) ≈ f_uw #oop, 3wrap + # test in-place version + unwrap!(f_wr, dims=1:3, circular_dims=(true, true, true)) + offset = first(f_uw) - first(f_wr) + @test (f_wr.+offset) ≈ f_uw #oop, 3wrap + end +end diff --git a/test/util.jl b/test/util.jl index b5976e243..7f5f14724 100644 --- a/test/util.jl +++ b/test/util.jl @@ -1,45 +1,5 @@ using DSP, Compat, Compat.Test -@testset "unwrap" begin - @test unwrap([0.1, 0.2, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] - @test unwrap([0.1, 0.2 + 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] - @test unwrap([0.1, 0.2 - 2pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] - @test unwrap([0.1, 0.2 - 2pi, 0.3 - 2pi, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] - @test unwrap([0.1 + 2pi, 0.2, 0.3, 0.4]) ≈ [0.1 + 2pi, 0.2 + 2pi, 0.3 + 2pi, 0.4 + 2pi] - @test unwrap([0.1, 0.2 + 6pi, 0.3, 0.4]) ≈ [0.1, 0.2, 0.3, 0.4] - - test_v = [0.1, 0.2, 0.3 + 2pi, 0.4] - res_v = unwrap(test_v) - @test test_v ≈ [0.1, 0.2, 0.3 + 2pi, 0.4] - res_v .= 0 - unwrap!(res_v, test_v) - @test res_v ≈ [0.1, 0.2, 0.3, 0.4] - @test test_v ≈ [0.1, 0.2, 0.3 + 2pi, 0.4] - unwrap!(test_v) - @test test_v ≈ [0.1, 0.2, 0.3, 0.4] - - # test multi-dimensional unwrapping - wrapped = [0.1, 0.2 + 2pi, 0.3, 0.4] - unwrapped = [0.1, 0.2, 0.3, 0.4] - wrapped = hcat(wrapped, wrapped) - unwrapped = hcat(unwrapped, unwrapped) - @test unwrap(wrapped, dims=2) ≈ wrapped - @test unwrap(wrapped, dims=1) ≈ unwrapped - @test unwrap!(copy(wrapped), dims=2) ≈ wrapped - @test unwrap!(copy(wrapped), dims=1) ≈ unwrapped - - # these should be implemented as part of #198 - @test_throws ArgumentError unwrap(wrapped, dims=1:2) - @test_throws ArgumentError unwrap!(similar(wrapped), wrapped, dims=1:2) - # this should eventually default to the above - @test_throws ArgumentError unwrap!(similar(wrapped), wrapped) - - # test unwrapping with other ranges - unwrapped = [1.0:100;] - wrapped = Float64[i % 10 for i in unwrapped] - @test unwrap(wrapped, range=10) ≈ unwrapped -end - @testset "hilbert" begin # Testing hilbert transform t = (0:1/256:2-1/256)