Hypercomplex Wavelets
From MMVLWiki
Contents |
Introduction
Complex wavelets are superior to real-valued wavelets because they are nearly shift-invariant. Complex wavelets yield amplitude-phase information in a similar way as the Fourier transform does. In contrast to the Fourier transform, wavelets allow to analyse the signal locally and thus can be applied to signals with a non-stationary statistic (such as images of a natural scene). In the same way as a one-dimensional signal requires complex numbers to represent the local structure of the signal, two-dimensional signals require hypercomplex numbers. Kingsbury has developed the Dual-Tree Hypercomplex Wavelet Transform (DHWT) which allows to recursively decompose a two-dimensional image.
Implementation
The implementation makes use of Selesnick's Hilbert transform pairs of wavelet bases. The implementation also requires the Ruby-extension HornetsEye which offers fast operations for n-dimensional arrays and hypercomplex numbers as element-types.
Under construction ...
#!/usr/bin/ruby
require 'selesnick'
require 'hornetseye'
require 'matrix'
require 'complex'
require 'hypercomplex'
include Hornetseye
class HWT
attr_reader(:h0)
attr_reader(:h0s)
attr_reader(:h1)
attr_reader(:h1s)
attr_reader(:g0)
attr_reader(:g0s)
attr_reader(:g1)
attr_reader(:g1s)
def initialize
# Change in text as well!
@h0, @h0s, @h1, @h1s, @g0, @g0s, @g1, @g1s = *hwt_wavelet_design( 5, 5, 4 )
@h0sx = MultiArray.to_multiarray( @h0s ).reshape( @h0s.size, 1 ).to_type( MultiArray::DFLOAT )
@g0sx = MultiArray.to_multiarray( @g0s ).reshape( @g0s.size, 1 ).to_type( MultiArray::DFLOAT )
@h1sx = MultiArray.to_multiarray( @h1s ).reshape( @h1s.size, 1 ).to_type( MultiArray::DFLOAT )
@g1sx = MultiArray.to_multiarray( @g1s ).reshape( @g1s.size, 1 ).to_type( MultiArray::DFLOAT )
@h0x = MultiArray.to_multiarray( @h0 ).reshape( @h0 .size, 1 ).to_type( MultiArray::DFLOAT )
@g0x = MultiArray.to_multiarray( @g0 ).reshape( @g0 .size, 1 ).to_type( MultiArray::DFLOAT )
@h1x = MultiArray.to_multiarray( @h1 ).reshape( @h1 .size, 1 ).to_type( MultiArray::DFLOAT )
@g1x = MultiArray.to_multiarray( @g1 ).reshape( @g1 .size, 1 ).to_type( MultiArray::DFLOAT )
@h0sy = MultiArray.to_multiarray( @h0s ).reshape( 1, @h0s.size ).to_type( MultiArray::DFLOAT )
@g0sy = MultiArray.to_multiarray( @g0s ).reshape( 1, @g0s.size ).to_type( MultiArray::DFLOAT )
@h1sy = MultiArray.to_multiarray( @h1s ).reshape( 1, @h1s.size ).to_type( MultiArray::DFLOAT )
@g1sy = MultiArray.to_multiarray( @g1s ).reshape( 1, @g1s.size ).to_type( MultiArray::DFLOAT )
@h0y = MultiArray.to_multiarray( @h0 ).reshape( 1, @h0 .size ).to_type( MultiArray::DFLOAT )
@g0y = MultiArray.to_multiarray( @g0 ).reshape( 1, @g0 .size ).to_type( MultiArray::DFLOAT )
@h1y = MultiArray.to_multiarray( @h1 ).reshape( 1, @h1 .size ).to_type( MultiArray::DFLOAT )
@g1y = MultiArray.to_multiarray( @g1 ).reshape( 1, @g1 .size ).to_type( MultiArray::DFLOAT )
end
def inspect
"HWT"
end
def downsamplex( n, c )
n.downsample( [2,1], [c,0] )
end
def upsamplex( n, c, size = -1 )
if size == -1
shape = nil
else
shape = [ size, n.shape[1] ]
end
n.upsample( [2,1], [c,0], shape )
end
def downsampley( n, c )
n.downsample( [1,2], [0,c] )
end
def upsampley( n, c, size = -1 )
if size == -1
shape = nil
else
shape = [ n.shape[0], size ]
end
n.upsample( [1,2], [0,c], shape )
end
def filtdx( img, filterx, c )
downsamplex( img.correlate( filterx ), c )
end
def filtdy( img, filtery, c )
downsampley( img.correlate( filtery ), c )
end
def filtux( img, filterx, c, size = -1 )
upsamplex( img, c, size ).correlate( filterx )
end
def filtuy( img, filtery, c, size = -1 )
upsampley( img, c, size ).correlate( filtery )
end
def hypercomplexarray( real, imag, jmag, kmag )
# real + imag * HyperComplex::I + ...
map = { MultiArray::SFLOAT => MultiArray::SHYPERCOMPLEX,
MultiArray::DFLOAT => MultiArray::DHYPERCOMPLEX }
map.default = MultiArray::SHYPERCOMPLEX
retval = MultiArray.new( map[ real.typecode ], *real.shape )
retval.real = real
retval.imag = imag
retval.jmag = jmag
retval.kmag = kmag
retval
end
def prepare( img )
# img * HyperComplex( 1.0, 1.0, 1.0, 1.0 )
hypercomplexarray( img, img, img, img )
end
def finalise( img )
0.25 * ( img.real + img.imag + img.jmag + img.kmag )
end
def decompose( wave, n = 1 )
c1=0
c2=0
real = wave.real
imag = wave.imag
jmag = wave.jmag
kmag = wave.kmag
realh0sx = filtdx( real, @h0sx, c1 )
imagg0sx = filtdx( imag, @g0sx, c1 )
jmagh0sx = filtdx( jmag, @h0sx, c1 )
kmagg0sx = filtdx( kmag, @g0sx, c1 )
realh1sx = filtdx( real, @h1sx, c1 )
imagg1sx = filtdx( imag, @g1sx, c1 )
jmagh1sx = filtdx( jmag, @h1sx, c1 )
kmagg1sx = filtdx( kmag, @g1sx, c1 )
real1 = filtdy( realh0sx, @h0sy, c1 )
imag1 = filtdy( imagg0sx, @h0sy, c1 )
jmag1 = filtdy( jmagh0sx, @g0sy, c1 )
kmag1 = filtdy( kmagg0sx, @g0sy, c1 )
real2 = filtdy( realh1sx, @h0sy, c1 )
imag2 = filtdy( imagg1sx, @h0sy, c1 )
jmag2 = filtdy( jmagh1sx, @g0sy, c1 )
kmag2 = filtdy( kmagg1sx, @g0sy, c1 )
real3 = filtdy( realh0sx, @h1sy, c2 )
imag3 = filtdy( imagg0sx, @h1sy, c2 )
jmag3 = filtdy( jmagh0sx, @g1sy, c2 )
kmag3 = filtdy( kmagg0sx, @g1sy, c2 )
real4 = filtdy( realh1sx, @h1sy, c2 )
imag4 = filtdy( imagg1sx, @h1sy, c2 )
jmag4 = filtdy( jmagh1sx, @g1sy, c2 )
kmag4 = filtdy( kmagg1sx, @g1sy, c2 )
img00 = hypercomplexarray( real1, imag1, jmag1, kmag1 )
img10 = hypercomplexarray( real2, imag2, jmag2, kmag2 )
img01 = hypercomplexarray( real3, imag3, jmag3, kmag3 )
img11 = hypercomplexarray( real4, imag4, jmag4, kmag4 )
if n == 1
[ img00, img10, img01, img11, wave.shape, [ c1, c2 ] ]
else
[ decompose( img00, n - 1 ), img10, img01, img11,
wave.shape, [ c1, c2 ] ]
end
end
def compose( wave, n = 1 )
if n == 1
img00 = wave[0]
else
img00 = compose( wave[0], n - 1 )
end
img10 = wave[1]
img01 = wave[2]
img11 = wave[3]
shape = wave[4]
c1 = wave[5][0]
c2 = wave[5][1]
real =
filtuy( filtux( img00.real, @h0x, c1, shape[0] ), @h0y, c1, shape[1] ) +
filtuy( filtux( img10.real, @h1x, c2, shape[0] ), @h0y, c1, shape[1] ) +
filtuy( filtux( img01.real, @h0x, c1, shape[0] ), @h1y, c2, shape[1] ) +
filtuy( filtux( img11.real, @h1x, c2, shape[0] ), @h1y, c2, shape[1] )
imag =
filtuy( filtux( img00.imag, @g0x, c1, shape[0] ), @h0y, c1, shape[1] ) +
filtuy( filtux( img10.imag, @g1x, c2, shape[0] ), @h0y, c1, shape[1] ) +
filtuy( filtux( img01.imag, @g0x, c1, shape[0] ), @h1y, c2, shape[1] ) +
filtuy( filtux( img11.imag, @g1x, c2, shape[0] ), @h1y, c2, shape[1] )
jmag =
filtuy( filtux( img00.jmag, @h0x, c1, shape[0] ), @g0y, c1, shape[1] ) +
filtuy( filtux( img10.jmag, @h1x, c2, shape[0] ), @g0y, c1, shape[1] ) +
filtuy( filtux( img01.jmag, @h0x, c1, shape[0] ), @g1y, c2, shape[1] ) +
filtuy( filtux( img11.jmag, @h1x, c2, shape[0] ), @g1y, c2, shape[1] )
kmag =
filtuy( filtux( img00.kmag, @g0x, c1, shape[0] ), @g0y, c1, shape[1] ) +
filtuy( filtux( img10.kmag, @g1x, c2, shape[0] ), @g0y, c1, shape[1] ) +
filtuy( filtux( img01.kmag, @g0x, c1, shape[0] ), @g1y, c2, shape[1] ) +
filtuy( filtux( img11.kmag, @g1x, c2, shape[0] ), @g1y, c2, shape[1] )
hypercomplexarray( real, imag, jmag, kmag )
end
end
Under construction ...
