Hypercomplex Wavelets
(Difference between revisions)
(Added introduction) |
(Added source-code) |
||
Line 5: | Line 5: | ||
way as a one-dimensional signal requires complex numbers to represent the local structure of the signal, two-dimensional signals | way as a one-dimensional signal requires complex numbers to represent the local structure of the signal, two-dimensional signals | ||
require hypercomplex numbers. | require hypercomplex numbers. | ||
+ | |||
+ | =Implementation= | ||
+ | [[Image:working.gif]] Under construction ... | ||
+ | <pre> | ||
+ | #!/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 | ||
+ | </pre> | ||
[[Image:working.gif]] Under construction ... | [[Image:working.gif]] Under construction ... |
Revision as of 14:55, 26 September 2007
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