,
Nick Kingsbury [ctb],
Rich Wareham [ctb]| Title: | Decimated and Undecimated 2D Complex Dual-Tree Wavelet Transform |
|---|---|
| Description: | An implementation of the decimated two-dimensional complex dual-tree wavelet transform as described in Kingsbury (1999) <doi:10.1098/rsta.1999.0447> and Selesnick et al. (2005) <doi:10.1109/MSP.2005.1550194>. Also includes the undecimated version and spectral bias correction described in Nelson et al. (2018) <doi:10.1007/s11222-017-9784-0>. The code is partly based on the 'dtcwt' Python library. |
| Authors: | Sebastian Buschow [aut, cre] ,
Nick Kingsbury [ctb],
Rich Wareham [ctb] |
| Maintainer: | Sebastian Buschow <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.5 |
| Built: | 2025-03-03 05:18:01 UTC |
| Source: | https://github.com/s6sebusc/dualtrees |
Matrices needed in order to eliminate the effect of "spectral leakage" for the local dtcwt-spectra.
A_b_bp A_bA_b_bp A_b
A list with entries N1024, N512, N256, N128, N64, N32, each containing the bias correction matrix of appropriate size.
As descibed in Nelson et al. (2018), the squared coefficients of the undecimated dtcwt (the local wavelet spectrum) can be used to obtain an estimate of local correlations in space. This estimator has a bias which mostly consists of an over-emphasis on the largest scales. It can be removed by multiplying each local spectrum with a matrix which depends only on the choice of wavelet and the dimensions of the field.
Calculated by brute force.
Nelson, J. D. B., A. J. Gibberd, C. Nafornita, and N. Kingsbury (2018) <doi:10.1007/s11222-017-9784-0>.
image( A_b_bp$N512 )image( A_b_bp$N512 )
A photograph of two meteorologists in front of the famous Japanese cherry trees in Bonn.
blossomblossom
A 256x256 matrix of gray-scale values
Armin Blanke, used with permission.
image(blossom, col=gray.colors(128,0,1))image(blossom, col=gray.colors(128,0,1))
Extend a matrix to the desired size.
pad(x, N, Ny = N, value = min(x, na.rm = TRUE)) put_in_mirror(x, N, Ny = N) period_bc(x, N, Ny = N)pad(x, N, Ny = N, value = min(x, na.rm = TRUE)) put_in_mirror(x, N, Ny = N) period_bc(x, N, Ny = N)
x |
a real matrix |
N |
the number of rows of the desired output |
Ny |
the number of columns of the desired output, defaults to N |
value |
the value with which the picture is padded by |
pad pads the fields with a constant value on all sides, be careful what you pick here. put_in_mirror reflects the input at all edges (with repeated end samples), period_bc simply repeats the input periodically. In any case, you can retrieve the initial area via bc$res[ bc$px, bc$py ].
a list containing the extended matrix ($res) and the positions of the original matrix within the extended one ($px and $py).
N and Ny must be at least as big as the input.
bc <- put_in_mirror( blossom, N=300 ) plot( bc ) print( range( bc$res[ bc$px, bc$py ] - blossom ) )bc <- put_in_mirror( blossom, N=300 ) plot( bc ) print( range( bc$res[ bc$px, bc$py ] - blossom ) )
Translate the centre of mass back and forth between polar and cartesian coordinates.
cen2xy(cen) xy2cen(xy)cen2xy(cen) xy2cen(xy)
cen |
the centre of mass of a wavelet spectrum (rho, phi, z), output of |
xy |
the centre of mass in cartesian coordinates (x, y, z), output of cen2xy |
These functions allow you to translate back and forth between the two coordinate systems. dt2cen represents the sepctrum's centre in cylinder coordinates because that is more intuitive than the x-y-z position within the hexagonal geometry. If you want to compare two spectra, it makes more sense to consider their distance in terms of x1-x2, y1-y2 since the difference in angle is only meaningful for reasonably large radii.
cen2xy is not the same thing as cen2uv !
transforms the angle and radius component of the spectral centre into vector components for visualization
cen2uv(cen)cen2uv(cen)
cen |
an |
The resulting vector field represents the degree of anisotropy and the dominant direction, averaged over all scales.
an nx x ny x 2 array, the two matrices representing the u- and v-component of the vector field.
dt <- fld2dt(blossom) ce <- dt2cen(dt) uv <- cen2uv(ce) uvplot( uv, z=blossom )dt <- fld2dt(blossom) ce <- dt2cen(dt) uv <- cen2uv(ce) uvplot( uv, z=blossom )
calculate the centre of mass of the local spectra in hexagonal geometry
dt2cen(dt, mask = NULL)dt2cen(dt, mask = NULL)
dt |
a |
mask |
a |
Each of the J x 6 spectral values is assigned a coordinate in 3D space with x(d,j)=cos(60*(d-1)), y(d,j)=sin(60*(d-1)), z(d,j)=j, where j denotes the scale and d the direction. Then the centre of mass in this space is calculated, the spectral values being the masses at each vertex. The x- and y-cooridnate are then transformed into a radius rho=sqrt(x^2+y^2) and an angle phi=15+0.5*atan2(y,x). rho measures the degree of anisotropy at each pixel, phi the orientation of edges in the image, and the third coordinate, z, the central scale. If a mask is provided, values where mask==TRUE are set to NA.
a nx x ny x 3 array where the third dimension denotes degree of anisotropy, angle and central scale, respectively.
Since the centre of mass is not defined for negative mass, any values below zero are removed at this point.
dt <- fld2dt(blossom) ce <- dt2cen(dt) image( ce[,,3], col=gray.colors(32, 0, 1) )dt <- fld2dt(blossom) ce <- dt2cen(dt) image( ce[,,3], col=gray.colors(32, 0, 1) )
average the output of fld2dt or dtcwt over space
dtmean(x)dtmean(x)
x |
either a J x nx x ny x 6 array of energies (output of dtcwt) or a list of complex wavelet coefficients (the output of |
a J x 6 matrix of spatially averaged energies
In the undecimated case, the coefficients are not averaged but summed up and then scaled by the area of the first level. This yields a comparable scale as the undecimated case.
These functions perform the dualtree complex wavelet analysis and synthesis, either with or without decimation.
dtcwt(fld, fb1 = near_sym_b, fb2 = qshift_b, J = NULL, dec = TRUE, verbose = FALSE) idtcwt(pyr, fb1 = near_sym_b, fb2 = qshift_b, verbose = TRUE)dtcwt(fld, fb1 = near_sym_b, fb2 = qshift_b, J = NULL, dec = TRUE, verbose = FALSE) idtcwt(pyr, fb1 = near_sym_b, fb2 = qshift_b, verbose = TRUE)
fld |
real matrix representing the field to be transformed |
fb1 |
A list of filter coefficients for the first level. Currently only |
fb2 |
A list of filter coefficients for all following levels. Currently only |
J |
number of levels for the decomposition. Defaults to |
dec |
whether or not the decimated transform is desired |
verbose |
if TRUE, the function tells you which level it is working on |
pyr |
a list containing arrays of complex coefficients for each level of the decomposition, produced by |
This is the 2D complex dualtree wavelet transform as described by Selesnick et al. (2005). It consists of four discerete wavelet transform trees, generated from two filter banks a and b by applying one set of filters to the rows and another (or the same) one to the columns. The 12 resulting coefficients are combined into six complex values representing six directions (15°, 45°, 75°, 105°, 135°, 165°).
In the decimated case (dec=TRUE), each convolution is followed by a downsampling by two, meaining that the size of the six coefficient fields is cut in half at each level. The decimated transform can be reversed to recover the original image. For the near_sym_b and qshift_b filter banks, this reconstrcution should be basically perfect. In the case of the the b_bp filters, non-negligible artifacts appear near +-45° edges.
if dec=TRUE a list of complex coefficient fields, otherwise a complex J x Nx x Ny x 6 array.
At present, the inverse transform only works if the input image had dimensions 2^N x 2^N. You can use boundaries to achieve that.
Nick Kingsbury (canonical MATLAB implementation), Rich Wareham (open source Python implementation, https://github.com/rjw57/dtcwt), Sebastian Buschow (R port).
Kingsbury, Nick (1999) <doi:10.1098/rsta.1999.0447>. Selesnick, I.W., R.G. Baraniuk, and N.C. Kingsbury (2005) <doi:10.1109/MSP.2005.1550194>
oldpar <- par( no.readonly=TRUE ) # forward transform dt <- dtcwt( blossom ) par( mfrow=c(2,3), mar=rep(2,4) ) for( j in 1:6 ){ image( blossom, col=grey.colors(32,0,1) ) contour( Mod( dt[[3]][ ,,j ] )**2, add=TRUE, col="green" ) } par( oldpar ) # exmaple for the inverse transform blossom_i <- idtcwt( dt ) image( blossom - blossom_i ) # example for a non-square case boy <- blossom[50:120, 50:150] bc <- put_in_mirror(boy, 128) dt <- dtcwt(bc$res) idt <- idtcwt(dt)[ bc$px, bc$py ]oldpar <- par( no.readonly=TRUE ) # forward transform dt <- dtcwt( blossom ) par( mfrow=c(2,3), mar=rep(2,4) ) for( j in 1:6 ){ image( blossom, col=grey.colors(32,0,1) ) contour( Mod( dt[[3]][ ,,j ] )**2, add=TRUE, col="green" ) } par( oldpar ) # exmaple for the inverse transform blossom_i <- idtcwt( dt ) image( blossom - blossom_i ) # example for a non-square case boy <- blossom[50:120, 50:150] bc <- put_in_mirror(boy, 128) dt <- dtcwt(bc$res) idt <- idtcwt(dt)[ bc$px, bc$py ]
Some of the filters implemented in the python package dtcwt.
qshift_b qshift_b_bp near_sym_b near_sym_b_bpqshift_b qshift_b_bp near_sym_b near_sym_b_bp
A list of high- and low-pass filters for analysis and synthesis
The near-sym filterbanks are biorthogonal wavelets used for the first level, they have 13 and 19 taps. The qshift filterbanks, each with 14 taps, are suitable for all higher levels of the dtcwt. The a- and b-filters form an approximate Hilbert-pair. The naming convention follows the python-package:
| h: | analysis |
| g: | synthesis |
| 0: | low-pass |
| 1: | high-pass |
| a,b: | shifted filters |
The b_bp-versions of the filterbanks contain a second high-pass for the diagonal directions, denoted by 2. They allow for better directional selectivity but prohibit perfect reconstruction.
'dtcwt' python package (https://github.com/rjw57/dtcwt)
Selesnick, I.W., R.G. Baraniuk, and N.C. Kingsbury (2005) <doi:10.1109/MSP.2005.1550194>
Kingsbury, N. (2006) <https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=7071567&isnumber=7065146>
Handles the transformation itself, boundary conditions and bias correction and returns the unbiased local wavelet spectrum at each grid-point.
fld2dt(fld, Nx = NULL, Ny = NULL, J = NULL, correct = TRUE, rsm = 0, verbose = FALSE, boundaries = "pad", fb1 = near_sym_b_bp, fb2 = qshift_b_bp)fld2dt(fld, Nx = NULL, Ny = NULL, J = NULL, correct = TRUE, rsm = 0, verbose = FALSE, boundaries = "pad", fb1 = near_sym_b_bp, fb2 = qshift_b_bp)
fld |
a real matrix |
Nx |
size to which the field is padded in x-direction |
Ny |
size to which the field is padded in y-direction |
J |
number of levels for the decomposition |
correct |
logical, whether or not to apply the bias correction |
rsm |
number of pixels to be linearly smoothed along each edge before applying the boundary conditions (see |
verbose |
whether or not you want the transform to talk to you |
boundaries |
how to handle the boundary conditions, either "pad", "mirror" or "periodic" |
fb1 |
filter bank for level 1 |
fb2 |
filter bank for all further levels |
The input is blown up to Nx x Ny and transformed by dtcwt(..., dec=FALSE). Then the original domain is cut out, the coefficients are squared and the bias is corrected (for details on the bias, see A).
an array of size J x nx x ny x 6 where dim(fld)=c(nx,ny)
Nelson, J. D. B., A. J. Gibberd, C. Nafornita, and N. Kingsbury (2018) <doi:10.1007/s11222-017-9784-0>
oldpar <- par( no.readonly=TRUE ) dt <- fld2dt( blossom ) par( mfrow=c(2,2), mar=rep(2,4) ) for( j in 1:4 ){ image( blossom, col=gray.colors(128, 0,1), xaxt="n", yaxt="n" ) for(d in 1:6) contour( dt[j,,,d], levels=quantile(dt[,,,], .995), col=d+1, add=TRUE, lwd=2, drawlabels=FALSE ) title( main=paste0("j=",j) ) } x0 <- seq( .1,.5,,6 ) y0 <- rep( 0.01,6 ) a <- .075 phi <- seq( 15,,30,6 )*pi/180 x1 <- x0 + a*cos( phi ) y1 <- y0 + a*sin( phi ) rect( min(x0,x1)-.05, min(y0,y1)-.05, max(x0,x1)+.05, max(y0,y1), col="black", border=NA ) arrows( x0, y0, x1, y1, length=.05, col=2:7, lwd=2, code=3 ) par( oldpar )oldpar <- par( no.readonly=TRUE ) dt <- fld2dt( blossom ) par( mfrow=c(2,2), mar=rep(2,4) ) for( j in 1:4 ){ image( blossom, col=gray.colors(128, 0,1), xaxt="n", yaxt="n" ) for(d in 1:6) contour( dt[j,,,d], levels=quantile(dt[,,,], .995), col=d+1, add=TRUE, lwd=2, drawlabels=FALSE ) title( main=paste0("j=",j) ) } x0 <- seq( .1,.5,,6 ) y0 <- rep( 0.01,6 ) a <- .075 phi <- seq( 15,,30,6 )*pi/180 x1 <- x0 + a*cos( phi ) y1 <- y0 + a*sin( phi ) rect( min(x0,x1)-.05, min(y0,y1)-.05, max(x0,x1)+.05, max(y0,y1), col="black", border=NA ) arrows( x0, y0, x1, y1, length=.05, col=2:7, lwd=2, code=3 ) par( oldpar )
square the wavelet coefficients and apply the bias correction
get_en(dt, correct = "none", N = ncol(dt))get_en(dt, correct = "none", N = ncol(dt))
dt |
array of |
correct |
type of correction, either |
N |
the smallest whole power of two larger than or equal to the dimensions of the input image, usually just |
The bias correction matrix should correspond to the filter bank used in the transform, for details on the matrices see A.
an array of the same dimensions as dt
Nelson, J. D. B., A. J. Gibberd, C. Nafornita, and N. Kingsbury (2018) <doi:10.1007/s11222-017-9784-0>
let a field decrease linearly towards its edges
smooth_borders(x, r)smooth_borders(x, r)
x |
a real matrix |
r |
a positive integer |
Values within the field are linearly reduced from their original value to the field minimum, starting r pixels away from the edge. This enforces truely periodic boundaries and removes sharp edges.
a matrix of the same dimensions as x
r must not be larger than min( dim(x) )/2.
image( smooth_borders(blossom, r=64), col=gray.colors(128,0,1) )image( smooth_borders(blossom, r=64), col=gray.colors(128,0,1) )
display the radial and angular component of the spectrum's centre as arrows.
uvplot(uv, z = NULL, x = NULL, y = NULL, col = "green", zcol = grDevices::gray.colors(32, 0, 1), n = 42, f = 1, length = 0.05, ...)uvplot(uv, z = NULL, x = NULL, y = NULL, col = "green", zcol = grDevices::gray.colors(32, 0, 1), n = 42, f = 1, length = 0.05, ...)
uv |
an array of dimension |
z |
image to show in the background, defaults to |
x, y
|
optional x- and y-coordinates for the plot, must match the dimensions of |
col |
color of the arrows |
zcol |
color scale for the image |
n |
number of arrows in one direction |
f |
factor by which to enlarge the arrows |
length |
length of the arrowhead in inches |
... |
further arguments passed to |
The pivot of the arrows is at the location to which the u- and v-component belong. No arrowhead is displayed since the egdges detcted by the cdtwt have an orientation but no sign. The default size of the arrows is such that a 'velocity' of 1 corresponds to 5% of the shorter image side.
uv <- cen2uv( dt2cen( fld2dt( blossom ) ) ) uvplot( uv, z=blossom )uv <- cen2uv( dt2cen( fld2dt( blossom ) ) ) uvplot( uv, z=blossom )