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.
Implementation
#!/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