fixed a problem with czt which causes wrap-around. This had an effect on resample_czt

Author: RainerHeintzmann

src/czt.jl

@@ -1,7 +1,7 @@
 export czt, iczt
 
 """
-    czt_1d(xin , scaled , d)
+    czt_1d(xin , scaled , d; kill_wrap=false, pad_value=zero(eltype(xin)))
 
 Chirp z transform along a single direction d of an ND array `xin` into the ND array 'xout'.
 Note that xin and xout can be the same array for inplace operations.
@@ -14,9 +14,11 @@ This code is based on a 2D Matlab version of the CZT, written by H. Gross et al.
 + `xin`: array to transform
 + `scaled`: factor to zoom into during the 1-dimensional czt. 
 + `d`: single dimension to transform (as a tuple)
++ `kill_wrap`: if true, the wrapped places will be set to pad_value
++ `pad_value`: the value that the wrap-around will be set to if `kill_wrap` is `true`. 
 
 """
-function czt_1d(xin, scaled, d)
+function czt_1d(xin, scaled, d; kill_wrap=false, pad_value=zero(eltype(xin)))
     sz=size(xin)
     # returns the real datatype
     rtype = real(eltype(xin))  
@@ -72,12 +74,17 @@ function czt_1d(xin, scaled, d)
             NDTools.slice(xout,d, o) .-= NDTools.slice(xin,d, 1) .* (1im).^mod(o-midp,4)
         end
     end
-    return xout
+    if kill_wrap && (scaled < 1.0)
+        nsz = Tuple(d == nd ? ceil(Int64, scaled * size(xin,d)) : size(xin,nd) for nd=1:ndims(xin))
+        return select_region(select_region(xout, new_size=nsz), new_size=size(xout), pad_value=pad_value)
+    else
+        return xout
+    end
     # xout .= g[dsize:(2*dsize-1)] .* reorient(fak, d, Val(ndims(xin)))
 end
 
 """
-    czt(xin , scale, dims=1:length(size(xin)))
+    czt(xin , scale, dims=1:length(size(xin)), kill_wrap=false)
 Chirp z transform of the ND array `xin`
 This code is based on a 2D Matlab version of the CZT, written by H. Gross.
 The tuple `scale` defines the zoom factors in the Fourier domain. Each has to be bigger than one.
@@ -90,6 +97,7 @@ The tuple `scale` defines the zoom factors in the Fourier domain. Each has to be
 + `xin`: array to transform
 + `scale`: a tuple of factors (one for each dimension) to zoom into during the czt. 
 + `dims`: a tuple of dimensions over which to apply the czt.
++ `kill_wrap`: if true, the wrapped places will be set to zero. Note that the `pad_value` argument is only allowed for 1d czts to not cause confusion.
 
 #Example:
 ```jdoctest
@@ -126,16 +134,16 @@ julia> zoomed = real.(ift(xft))
   0.0239759   -0.028264    0.0541186  -0.0116475   -0.261294   0.312719  -0.261294  -0.0116475    0.0541186  -0.028264
 ```
 """
-function czt(xin::Array{T,N}, scale, dims=1:length(size(xin)))::Array{complex(T),N} where {T,N}
+function czt(xin::Array{T,N}, scale, dims=1:length(size(xin)); kill_wrap=false)::Array{complex(T),N} where {T,N}
     xout = xin
     for d in dims
-        xout = czt_1d(xout, scale[d], d)
+        xout = czt_1d(xout, scale[d], d; kill_wrap=kill_wrap)
     end
     return xout
 end
 
 """
-    iczt(xin , scale, dims=1:length(size(xin)))
+    iczt(xin , scale, dims=1:length(size(xin)); kill_wrap=false)
 Inverse chirp z transform of the ND array `xin`
 This code is based on a 2D Matlab version of the CZT, written by H. Gross.
 The tuple `scale` defines the zoom factors in the Fourier domain. Each has to be bigger than one.
@@ -146,6 +154,7 @@ The tuple `scale` defines the zoom factors in the Fourier domain. Each has to be
 + `xin`: array to transform
 + `scale`: a tuple of factors (one for each dimension) of the the inverse czt. 
 + `dims`: a tuple of dimensions over which to apply the inverse czt.
++ `kill_wrap`: if true, the wrapped places will be set to zero. Note that the `pad_value` argument is only allowed for 1d czts to not cause confusion.
 
 #See also: czt, czt_1d
 
@@ -184,6 +193,6 @@ julia> iczt(xft,(1.2,1.3))
  -0.0965531+0.0404296im  -0.159713+0.0637132im    0.48095+0.0775406im    0.67753-0.263814im        0.77553-0.121603im    0.660335-0.00736904im   0.495205-0.135059im   -0.163859+0.125535im
  ```
 """
-function iczt( xin , scale, dims=1:length(size(xin)))
-    conj(czt(conj(xin), scale, dims)) / prod(size(xin))
+function iczt( xin , scale, dims=1:length(size(xin)); kill_wrap=false)
+    conj(czt(conj(xin), scale, dims; kill_wrap=kill_wrap)) / prod(size(xin))
 end

src/resampling.jl

@@ -167,16 +167,27 @@ function upsample2_abs2(mat::AbstractArray{T, N}; dims=1:N) where {T,N}
 end
 
 """
-    resample_czt(arr, rel_zoom; shear=nothing, shear_dim=nothing, fix_nyquist=false, new_size = size(arr), rel_pad=0.2)
+    resample_czt(arr, rel_zoom; shear=nothing, shear_dim=nothing, fix_nyquist=false, new_size = size(arr), do_damp=false, rel_pad=0.2, kill_wrap=true)
 
 resamples the image with fixed factors or a list of separable functions using the chirp z transform algorithm.
 The data is first padded by a relative amount `rel_pad` which is needed to avoid wrap-around problems.
 As opposed to `resample()`, this routine allows for arbitrary non-integer zoom factors.
-It is reasonably fast but only allows a stretch (via `rel_zoom`) and a shift (via `shear` in pixels) per line or column
+It is reasonably fast but only allows a stretch (via `rel_zoom`) and a shift (via `shear` in pixels) per line or column.
 
 Note that each entry of the tuple in `rel_zoom` or `shear` describes the zoom or shear to apply to all other dimensions individually
 per entry along this dimension number. 
 
+# Arguments:
++ `arr`: array to resample
++ `rel_zoom`: factors to zoom as a tuple or a tuple of functions defining the zooms
++ `shear`: a tuple of shears or a tuple of shear functions defining the shears
++ `shear_dim`: which dimension to shear
++ `fix_nyquist`: defines whether to apply `fix_nyquist` when using the apply_shift_strength! function.
++ `do_damp`: applies a padding and damping outside the region  to zoom, to avoid artefacts
++ `rel_pad`: amount of padding to apply, if `do_damp` is true
++ `kill_wrap`: removes the wrap-around when zooming out.
++ `new_size`: size of the result array. If not provided the same as the input size will be used.
+
 # Examples
 ```jdoctest
 julia> using TestImages, NDTools, View5D
@@ -192,21 +203,28 @@ julia> d = resample_czt(a, (x->1.0,x->1.0), shear=(x->50*x^2,0.0)); # a more com
 julia> @ve a,b,c,d # visualize distortions
 ```
 """
-function resample_czt(arr::AbstractArray{T,N}, rel_zoom; shear=nothing, shear_dim=nothing, fix_nyquist=false, new_size = size(arr), rel_pad=0.2, do_damp=false, center=CtrMid) where {T,N}
+function resample_czt(arr::AbstractArray{T,N}, rel_zoom; shear=nothing, shear_dim=nothing, fix_nyquist=false, new_size = size(arr), rel_pad=0.2, do_damp=false, center=CtrMid, kill_wrap=true, pad_value=zero(eltype(arr))) where {T,N}
     RT = real(T)
     orig_size = size(arr)
     if do_damp
         arr = damp_edge_outside(arr, rel_pad)
     else
         arr = copy(arr)
     end
-    for d in 1:length(rel_zoom)
+    for d in eachindex(rel_zoom)
         sd = mod(d,ndims(arr))+1
         if !isnothing(shear_dim)
             sd = shear_dim[d]
         end
         my_zoom = 1.0
-        # case of a list of zoom numbers
+        # resize the array before zooming, in case the result array is larger
+        # This may be a bit wasteful depending on the zoom factors, but this case also
+        # covers the zoomed shearing case
+        if (new_size[d] > size(arr,d))
+            nz = Tuple(d == nd ? new_size[nd] : size(arr,nd) for nd=1:ndims(arr))
+            arr = collect(select_region(arr, new_size=nz))
+        end
+    # case of a list (or tuple) of zoom numbers
         if (isa(rel_zoom, Tuple) && isa(rel_zoom[1], Number)) || isa(rel_zoom, Number)
             my_zoom = rel_zoom[d]
             f_res = ift(arr, d)
@@ -219,9 +237,9 @@ function resample_czt(arr::AbstractArray{T,N}, rel_zoom; shear=nothing, shear_di
                 FourierTools.apply_shift_strength!(f_res, f_res, shifts, d, sd, -myshear, fix_nyquist)
             end
             if T<:Real
-                f_res = real(FourierTools.czt_1d(f_res, my_zoom, d))
+                f_res = real(FourierTools.czt_1d(f_res, my_zoom, d; kill_wrap=kill_wrap, pad_value=pad_value))
             else
-                f_res = T.(FourierTools.czt_1d(f_res, my_zoom, d))
+                f_res = T.(FourierTools.czt_1d(f_res, my_zoom, d; kill_wrap=kill_wrap, pad_value=pad_value))
             end
             select_region!(f_res, arr) 
         # case of position dependent zoom functions
@@ -251,9 +269,9 @@ function resample_czt(arr::AbstractArray{T,N}, rel_zoom; shear=nothing, shear_di
                 end
                 f_res = let 
                     if T<:Real
-                        real(FourierTools.czt_1d(f_res, my_zoom, 1))
+                        real(FourierTools.czt_1d(f_res, my_zoom, 1; kill_wrap=kill_wrap, pad_value=pad_value))
                     else
-                        FourierTools.czt_1d(f_res, my_zoom, 1)
+                        FourierTools.czt_1d(f_res, my_zoom, 1; kill_wrap=kill_wrap, pad_value=pad_value)
                     end
                 end
                 select_region!(f_res, slice)

test/czt.jl

@@ -10,6 +10,13 @@ using NDTools # this is needed for the select_region! function below.
         @test ≈(iczt(czt(y,zoom),zoom),  y, rtol=1e-5)
         zoom = (2.0,2.0)
         @test ≈(czt(y,zoom),  select_region(upsample2(ft(y), fix_center=true),new_size=size(y)), rtol=1e-5)
+        # zoom smaller 1.0 causes wrap around:
+        zoom = (0.5,2.0)
+        @test abs(czt(y,zoom)[1,1]) > 1e-5
+        zoom = (0.5,2.0)
+        # check if the kill_wrap works
+        @test abs(czt(y,zoom; kill_wrap=true)[1,1]) == 0.0
+        @test abs(iczt(y,zoom; kill_wrap=true)[1,1]) == 0.0
         # @vt czt(y,zoom)  select_region(upsample2(ft(y), fix_center=true), new_size=size(y))
     end
 end