Resampling

resample(arr, new_size [, normalize=true])

Calculates the sinc interpolation of an arr on a new array size new_size. It is a re-evaluation of the Fourier series at new grid points. new_size can be arbitrary. Each dimension is then independently either up or downsampled.

This method is based on FFTs and therefore implicitly assumes periodic boundaries and a finite frequency support.

normalize=true by default multiplies by an appropriate factor so that the array size is included in the scaling. This results in an array having roughly the same mean intensity.

Basic Principle

If size(new_size)[i] > size(arr)[i], we apply zero padding in Fourier space.

If size(new_size)[i] < size(arr)[i], we cut out a centered part of the Fourier spectrum.

We apply some tricks at the boundary to increase accuracy of highest frequencies.

For real arrays we use rfft based operations, for complex one we use fft based ones.

Examples

sinc interpolation of 2 datapoints result in an approximation of cosine.

julia> resample([2.0, 0.0], (6,))
6-element Vector{Float64}:
 2.0
 1.5
 0.5
 0.0
 0.5
 1.5

julia> resample([2.0, 0.0], (6,)) ≈ 1 .+ cos.(2π .* (0:5)./6)
true

resample_by_FFT(mat, new_size)

Does a resampling based on fft. This function is called by resampling.

resample_by_RFFT(mat, new_size)

Does a resampling based on rfft. This function is called by resampling.

resample_nfft(img, new_pos, dst_size=nothing; pixel_coords=false, is_local_shift=false, is_src_coords=true, reltol=1e-9)

resamples an ND-array to a set of new positions new_pos measured in either in pixels (pixel_coords=true) or relative (Fourier-) image coordinates (pixel_coords=false). new_pos can be

• an array of Tuples specifying the zoom along each direction
• an N+1 dimensional array (for N-dimensional input data img) of destination positions, the last dimension enumerating the respective destination coordinate dimension.
• a function accepting a coordinate Tuple and yielding a destination position Tuple.

resample_nfft can perform a large range of possible resamplings. Note that the default setting is is_src_coords=true which means that the source coordinates of each destination position have to be specified. This has the advantage that the result has usually less artefacts, but the positions may be more less convenient to specify.

Arguments

• img: the image to apply resampling to
• new_pos`: specifies the resampling. See description above.
• dst_size: this argument optionally defines the output size. If you require a different result size for new_pos being a function or with is_src_coords=true, state it here. By default (dst_size=nothing) the destination size will be inferred form the argument new_pos or assumed to be size(img).
• is_local_shift: specifies, whether the resampling coordinates refer to a relative shift or absoluter coordinates
• is_in_pixels: specifies whether the coordinates (or relative distances) are given in pixel pitch units (is_in_pixels=true) or in units relative to the array sizes (Fourier convention)
• is_src_coords: specifies, whether the resampling positions refer to sampling at source (is_src_coords=true) or destination coordinates
• reltol: will be used as an argument to the nfft function spedifying the relative precision to calculate to

See also: resample, resample_czt

Examples

julia> using FourierTools, TestImages, NDTools, View5D, IndexFunArrays

julia> a = Float32.(testimage("resolution"));

julia> b = resample_nfft(a, t -> (2.5f0 *sign(t[1])*t[1]^2, t[2]*(0.5f0+t[1]))); # a complicated deformation

julia> sz = size(a);

# stacking only the displacement along the last dimension:
julia> new_pos = cat(xx(sz,scale=ScaFT), zeros(sz), dims=3);

julia> c = resample_nfft(a, new_pos, is_local_shift=true); # stretch along x using an array

julia> new_pos = cat(.-xx(sz,scale=ScaFT)./2, zeros(sz), dims=3);

julia> c2 = resample_nfft(a, new_pos, is_local_shift=true, is_src_coords=false); # stretch along x using an array

# Notice the difference in brightness between c and c2
julia> @ve a b c c2 # visualize distortion and x-shrinks. 

# Lets try a 2D rotation:
# define a rotation operation
julia> rot_alpha(a, t) = (cosd(a)*t[1] + sind(a)*t[2], -sind(a)*t[1]+cosd(a)*t[2])

# positions as an array of tuples
julia> new_pos = rot_alpha.(10.0, idx(a, scale=ScaFT))

# lets do the resampling, this time by specifying the destination coordinates:
julia> d = resample_nfft(a, new_pos, is_src_coords=false);

#display the result
julia> @ve a d

#how about a spiral deformation?
julia> new_pos = rot_alpha.(rr(a), idx(a, scale=ScaFT))

julia> e = resample_nfft(a, new_pos);

julia> f = resample_nfft(a, new_pos, is_src_coords=false);

# observe the artefacts generated by undersampling in the destination grid
julia> @ve a e f

resample_czt(arr, rel_zoom; shear=nothing, shear_dim=nothing, fix_nyquist=false, new_size = size(arr), 
             do_damp=false, rel_pad=0.2, remove_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.

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 applyshiftstrength! 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
• remove_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

julia> using TestImages, NDTools, View5D

julia> a = Float32.(testimage("resolution"));

julia> b = resample_czt(a, (0.99,1.3)); # just zooming and shrinking

julia> c = resample_czt(a, (x->0.8+x^2/3,x->1.0+0.2*x)); # a more complicated distortion

julia> d = resample_czt(a, (x->1.0,x->1.0), shear=(x->50*x^2,0.0)); # a more complicated distortion

julia> @ve a,b,c,d # visualize distortions

upsample2(mat::AbstractArray{T, N}; dims=1:N, fix_center=false, keep_singleton=false)

Upsamples by a factor of two in all dimensions. The code is optimized for speed by using subpixelshifts rather than Fourier resizing. By default the first pixel maintains its position. However, this leads to a shift of the center (size[d]÷2+1) in the resulting array for uneven array sizes. fix_center=true can be used to remedy this and the result array center position will agree to the source array center position. keep_singleton=true will not upsample dimensions of size one. Note that upsample2 is based on Fourier-shifting and you may have to deal with wrap-around problems.

julia> upsample2(collect(collect(1.0:9.0)'))
2×18 Matrix{Float64}:
 1.0  0.24123  2.0  3.24123  3.0  2.93582  4.0  5.0  5.0  5.0  6.0  7.06418  7.0  6.75877  8.0  9.75877  9.0  5.0
 1.0  0.24123  2.0  3.24123  3.0  2.93582  4.0  5.0  5.0  5.0  6.0  7.06418  7.0  6.75877  8.0  9.75877  9.0  5.0

julia> upsample2(collect(collect(1.0:9.0)'); fix_center=true, keep_singleton=true)
1×18 Matrix{Float64}:
 5.0  1.0  0.24123  2.0  3.24123  3.0  2.93582  4.0  5.0  5.0  5.0  6.0  7.06418  7.0  6.75877  8.0  9.75877  9.0

upsample2_abs2(mat::AbstractArray{T, N}; dims=1:N)

Upsamples by a factor of two and applies the abs2 operation. The code is optimized for speed.

barrel_pin(arr, rel=0.5)

emulates a barrel (rel>0) or a pincushion (rel<0) distortion. The distortions are calculated using resample_czt() with separable quadratic zooms.

See also: resample_czt()

Examples

julia> using TestImages, NDTools, View5D

julia> a = Float32.(testimage("resolution"))

julia> b = barrel_pin(a,0.5) # strong barrel distortion

julia> c = barrel_pin(a,-0.5) # strong pin-cushion distortion

julia> @ve a,b,c # visualize distortions