Two-Dimensional Fourier Transform And Linear Filtering
Two-Dimensional Fourier Transformand Linear FilteringYao WangPolytechnic School of Engineering, New York University Yao Wang, 2016EL-GY 6123: Image and Video Processing1
Outline General concept of signals and transforms– Representation using basis functions Continuous Space Fourier Transform (CSFT)– 1D - 2D– Concept of spatial frequency Discrete Space Fourier Transform (DSFT) and DFT– 1D - 2D Continuous space convolution Discrete space convolution Convolution theorem Yao Wang, 2016EL-GY 6123: Image and Video Processing2
Signal in 2D Space General 2D continuous space signal: f(x,y)– Can have infinite support: x,y (-infty, , infty)– f(x,y) can generally take on complex values General 2D discrete space signal: f(m,n)– Can have infinite support: m,n -infty, , 0, 1,., infty– f(m,n) can generally take on complex values Each color component of an image is a 2D real signal with finitesupport–––– Yao Wang, 2016MxN image: m 0,1, ,M-1, n 0,1, ,N-1We will use first index for row, second index for columnWe will consider a single color component onlySame operations can be applied to each componentEL-GY 6123: Image and Video Processing3
Transform Representation of Signals Transforms are decompositions of a function f(x) into some basisfunctions Ø(x, u). u is typically the freq. index. Yao Wang, 2016EL-GY 6123: Image and Video Processing4
Illustration of Decomposition in Vector SpaceΦ3fα3f α1Φ1 α2Φ2 α3Φ3oα1α2Φ2Φ1 Yao Wang, 2016EL-GY 6123: Image and Video Processing5
Decomposition of 1D Signal Ortho-normal basis function 1, u1 u2 φ ( x, u1 )φ * ( x, u2 )dx 0, u1 u2 Forward transformF (u ) f ( x),φ ( x, u ) f ( x)φ * ( x, u )dxProjection off(x) onto ϕ(x,u) Inverse transformf ( x) F (u )φ ( x, u )duRepresenting f(x) as sum of φ(x,u)for all u, with weight F(u) Yao Wang, 2016EL-GY 6123: Image and Video Processing6
1D Continuous Time Fourier Transform Basis functionφ ( x, u) ej 2πux, u ( , ). Forward Transform F (u ) F{ f ( x)} f ( x)e j 2πux dx Inverse Transform f ( x) F {F (u )} F (u )e j 2πux du 1 Yao Wang, 2016EL-GY 6123: Image and Video Processing7
Important Transform Pairsf ( x) 1 f ( x) e j 2πf 0 xF (u ) δ (u )F (u ) δ (u f 0 )1f ( x) cos( 2πf 0 x) F (u ) (δ (u f 0 ) δ (u f 0 ) )21(δ (u f 0 ) δ (u f 0 ) )f ( x) sin(2πf 0 x) F (u ) 2j 1, x x0f ( x) 0, otherwise Yao Wang, 2016F (u ) sin(2πx0u ) 2 x0 sinc(2 x0u )πusin(πt )where, sinc(t ) πtEL-GY 6123: Image and Video Processing8
FT of the Rectangle Functionsin(2πx0u )sin(πt )F (u ) 2 x0 sinc( 2 x0u ) where , sinc( t ) πuπtf(x)-1f(x)x0 11x-2x0 22xNote first zero occurs at u0 1/(2 x0) 1/pulse-width, other zeros are multiples of this. Yao Wang, 2016EL-GY 6123: Image and Video Processing9
IFT of Ideal Low Pass Signal What is f(x)?F(u)-u0 Yao Wang, 2016u0uEL-GY 6123: Image and Video Processing10
Representation of FT Generally, both f(x) and F(u) are complex Two representations– Real and Imaginary– Magnitude and PhaseF (u ) R(u ) jI (u )F (u ) A(u )e jφ (u ) , whereI (u )A(u ) R(u ) I (u ) , φ (u ) tanR(u )22II(u)F(u) 1Φ(u)RR(u) RelationshipR(u ) A(u ) cos φ (u ), I (u ) A(u ) sin φ (u ) Power spectrumP(u) A(u) F (u) F (u) F (u)2 Yao Wang, 2016EL-GY 6123: Image and Video Processing*211
What if f(x) is real? Real world signals f(x) are usually realF(u) is still complex, but has special propertiesF * (u ) F ( u )R(u ) R( u ), A(u ) A( u ), P(u ) P( u ) : even functionI (u ) I ( u ),φ (u ) φ ( u ) : odd function Yao Wang, 2016EL-GY 6123: Image and Video Processing12
Property of Fourier Transform Duality Linearity Scalingf (t ) F (u )F (t ) f ( u )F{a1 f1 ( x) a2 f 2 ( x)} a1F{ f1 ( x)} a2 F{ f 2 ( x)}F{af ( x)} aF{ f ( x)} Translationf ( x x0 ) F (u )e j 2πx0u , Convolutionf ( x)e j 2πu0 x F (u u0 )f ( x) g ( x) f ( x α ) g (α )dαf ( x) g ( x) F (u )G(u )We will review convolution later! Yao Wang, 2016EL-GY 6123: Image and Video Processing13
Two Dimension Continuous Space FourierTransform (CSFT) Basis functionsφ ( x, y; u, v) e j (2πux 2πvy ) e j 2πuxe j 2πvy , u, v ( , ). Forward – TransformF (u, v) F{ f ( x, y)} f ( x, y)e j 2π (ux vy ) dxdy Inverse – Transformf ( x, y ) F {F (u, v)} 1 F (u, v)e j 2π (ux vy ) dudv– Representing a 2D signal as sum of 2D complex exponentialsignals Yao Wang, 2016EL-GY 6123: Image and Video Processing14
Illustration of Spatial Frequencyf (x , y) sin(10π x)f x 5, f y 0, f m 5,φ 0f x : !unit !1!cycle/image7width!f y : unit 1!cycle/image7heightf (x , y) sin(10π x 20π y)f x 5, f y 10, f m 125,φ atan( 2)!! Yao Wang, 2016EL-GY 6123: Image and Video Processing15
Example 1f ( x, y ) sin 4πx cos 6πyF{sin 4πx} sin 4πxe j 2π (ux vy ) dxdyf(x,y) sin 4πxe j 2πux dx e j 2πvy dy sin 4πxe j 2πux dxδ (v)1 (δ (u 2) δ (u 2))δ (v)2j1 (δ (u 2, v) δ (u 2, v))2j , x y 0where δ ( x, y ) δ ( x)δ ( y ) 0, otherwiseuF(u,v)v1Likewise, F{cos 6πy} (δ (u , v 3) δ (u, v 3))2 Yao Wang, 2016EL-GY 6123: Image and Video Processing16
Example 2f ( x, y ) sin( 2πx 3πy ) {}(1 j ( 2πx 3πy ) j ( 2πx 3πy )e e2j)F e j ( 2πx 3πy ) e j ( 2πx 3πy ) e j 2π ( xu yv ) dxdy e j 2πx e j 2πux dx e j 3πy e j 2πyv dy33 δ (u 1)δ (v ) δ (u 1, v )223Likewise, F e j ( 2πx 3πy ) δ (u 1, v )2Therefore,{F {sin( 2πx 3πy )} Yao Wang, 2016}1 33 δ (u 1, v ) δ (u 1, v ) 2j 22 EL-GY 6123: Image and Video Processing[X,Y] meshgrid(-2:1/16:2,-2:1/16:2);f sin(2*pi*X 3*pi*Y);imagesc(f); colormap(gray)Truesize, axis off;uF(u,v)v17
Properties of 2D FT (1) LinearityF{a1 f1 ( x, y) a2 f 2 ( x, y)} a1F{ f1 ( x, y)} a2 F{ f 2 ( x, y)} Translationf ( x x0 , y y0 ) F (u, v)e j 2π ( x0u y0v ) ,f ( x, y)e j 2π (u0 x v0 y ) F (u u0 , v v0 ) Conjugationf * ( x, y) F * ( u, v) Yao Wang, 2016EL-GY 6123: Image and Video Processing18
Properties of 2D FT (2) Symmetryf ( x, y) is real F (u, v) F ( u, v) Convolution– Definition of convolutionf ( x, y) g ( x, y) f ( x α , y β ) g (α , β )dαdβ– Convolution theoryf ( x, y) g ( x, y) F (u, v)G(u, v)We will describe 2D convolution later! Yao Wang, 2016EL-GY 6123: Image and Video Processing19
Separability of 2D FT and Separable Signal Separability of 2D FTF2 { f ( x, y )} Fy {Fx { f ( x, y )}} Fx {Fy { f ( x, y )}}– where Fx, Fy are 1D FT along x and y.– one can do 1DFT for each row of original image,then 1D FT along each column of resulting image Separable Signal– f(x,y) fx(x)fy(y)– F(u,v) Fx(u)Fy(v), where Fx(u) Fx{fx(x)}, Fy(u) Fy{fy(y)}– For separable signal, one can simply compute two1D transforms and take their product! Yao Wang, 2016EL-GY 6123: Image and Video Processing20
Example 1f ( x, y ) sin(3πx) cos(5πy )f x ( x) sin(3πx) f y ( y ) cos(5πy ) 1(δ (u 3 / 2) δ (u 3 / 2))2j1Fy (v) (δ (v 5 / 2) δ (v 5 / 2) )2Fx (u ) F (u, v) Fx(u ) Fy(v) 1 35353535 δ (u , v ) δ (u , v ) δ (u , v ) δ (u , v ) 4j 22222222 uv Yao Wang, 2016EL-GY 6123: Image and Video Processing21
Example 2 1, x x0 , y y0f ( x, y ) 0, otherwiseF (u, v) 4 x0 y0 sinc(2 x0u ) sinc(2 y0v)x2ux0 2y0 1-11vy-2 Yao Wang, 2016w/EL-GY 6123: Image and Video Processinglogrithmic mapping22
Example 3: Gaussian Signal Still a Gaussian Function! 1D Gaussian Signal 2D Gaussian Signal!! !!!!!!exp exp ! exp !2! !2!2!!!!!!! ! !1 exp ! exp ! exp ,! !!2!2!2!2!"! Note that STD β in freq. inversely related to STD σ inspace Yao Wang, 2016EL-GY 6123: Image and Video Processing23
Illustration of Gaussian Signalσ 1β 0.16g(x)G(u)σ 2 Yao Wang, 2016β 0.08EL-GY 6123: Image and Video Processing24
2D Guaussian FunctionPlot as imageSurface plot (surf( )) Yao Wang, 2016EL-GY 6123: Image and Video Processing25
Important Transform Pairs All following signals are separable and can be provedby applying the separability of the CSFTf (x, y) 1 F(u, v) δ (u, v), whereδ (u, v) δ (u)δ (v)f (x) e j 2 π ( f1x f2 y) F(u) δ (u f1, v f2 )1f (x, y) cos(2π ( f1 x f2 y)) F(u) (δ (u f1 , v f2 ) δ (u f1 , v f2 ))21f (x, y) sin(2π ( f1 x f2 y)) F(u) (δ (u f1 , v f2 ) δ (u f1 , v f2 ))2j# 1, x x , y y%00f (x, y) %& 0, otherwisesin(2π x0 u) sin(2π y0 v)F(u, v) 4x0 y0 sinc(2x0 u)sinc(2y0 v)πuπvsin(π t)where sinc(t) πt Yao Wang, 2016EL-GY 6123: Image and Video Processing Constant ó impulse at(0,0) freq. Complexexponentialó impulse ata particular2D freq. 2D box ó 2Dsinc function Gaussian ó Gaussian26
Rotation Let x r cosθ , y r sin θ , u ρ cos ω , 2D FT in polar coordinate (r, θ) and (ρ, Ø)F (ρ ,φ ) 0 2π0f ( r , θ )ev ρ sin ω. j 2π ( r cos θ ρ cos φ r sin θ ρ sin φ )rdrdθ f (r , θ )e j 2πrρ cos(θ φ ) rdrdθ Propertyf (r , θ θ 0 ) F ( ρ , φ θ 0 ) Proof: Homework! Yao Wang, 2016EL-GY 6123: Image and Video Processing27
Example of Rotation Yao Wang, 2016EL-GY 6123: Image and Video Processing28
1D Fourier Transform For Discrete TimeSequence (DTFT) (Review) f(n) is a 1D discrete time sequenceForward Transform F (u ) f (n)e j 2πun1/ 2F (u )e j 2πun dun Inverse Transformf ( n) 1/ 2 Representing f(n) as weighted sum of many complex sinusoidalsignals with frequency u, F(u) is the weight F(u) indicate the “amount” of sinusoidal component with freq. u insignal f(n) u digital frequency the number of cycles per integer sample.Period 1/ u (must be equal or greater than 2 samples- u 1/2) Yao Wang, 2016EL-GY 6123: Image and Video Processing29
Properties unique for DTFT Periodicity– F(u) F(u 1)– The FT of a discrete time sequence is only considered for u є (½ , ½), and u ½ is the highest discrete frequency Symmetry for real sequencesf (n) f * (n) F (u ) F * ( u ) Yao Wang, 2016 F (u ) F ( u ) F (u ) is symmetricEL-GY 6123: Image and Video Processing30
Example 1, n 0,1,., N 1;f ( n) 0, others F (u ) N 1 e j 2πnu 1 e j 2πuN1 en 0 j 2πu e jπ ( N 1)usin 2πu ( N / 2)sin 2πu (1 / 2)f(n)0N-1N 10nThere are N/2 zeros in (0, ½], 1/N apart Yao Wang, 2016EL-GY 6123: Image and Video Processing31
1D Discrete Fourier Transform (DFT) DSFT of N-pt sequenceN 1 N-pt DFT of N-pt sequenceN 1F(u) f (n)e j2π unn 0!FN (k) f (n)en 0! j2πknN k F u N DFT is the sampled version of DTFT with samples at 0,1/N, , (N-1)/N. FFT: Fast algorithm for computing DFT– Direct computation takes N 2 operations– FFT takes N log(N) operations! Yao Wang, 2016EL-GY 6123: Image and Video Processing32
Discrete Space Fourier Transform (DSFT) forTwo Dimensional Signals Let f(m,n) represent a 2D sequence Forward TransformF (u, v) f (m, n)e j 2π ( mu nv )m n Inverse Transform1/ 2f (m, n) 1/ 2 1/ 2 1/ 2F (u, v)e j 2π ( mu nv ) dudv Representing an image as weighted sum of many 2D complexsinusoidal images u: number of cycles per vertical sample (vertical freq.) v: number of cycles per horizontal sample (horizontal freq) Yao Wang, 2016EL-GY 6123: Image and Video Processing33
Spatial Frequency for Digital Imagesf s 125, ϕ arctan(2)If the image has 256x256 pixels, fx 5 cycles per width (analog frequency) - u 5cycles/256 pixels 5/256 cycles/sample (digital frequency)When both horizontal and vertical frequency are non-zero, we see directionalpatterns. fs is the frequency along the direction with the maximum change(orthogonal to the lines)Note that 1/u may not correspond to integer.If u a/b, digital period b. Eg u 3/8, Digital Period 8 Yao Wang, 2016EL-GY 6123: Image and Video Processing34
PeriodicityF (u, v) f (m, n)e j 2π ( mu nv )m n F(u,v) is periodic in u, v with period 1, i.e., for allintegers k, l:– F(u k, v l) F(u, v) To see this considere Yao Wang, 2016 j 2π ( m ( u k ) n ( v l )) e j 2π ( mu nv ) j 2π ( mk nl ) e j 2π ( mu nv )eEL-GY 6123: Image and Video Processing35
Example: Delta Function Fourier transform of a delta functionF (u, v) δ (m, n)e j 2π ( mu nv) 1m n δ (m, n) 1– Delta function contains all frequency components with equalweights! Inverse Fourier transform of a delta functionf (m, n) 0.5 0.5 0.5 0.51 δ (u, v) Yao Wang, 2016δ (u, v)e j 2π ( mu nv ) dudv 1EL-GY 6123: Image and Video Processing36
Example21 1f (m, n) 000 1 2 1 mf(m,n)nF (u, v) 1e j 2π ( 1*u 1*v ) 2e j 2π ( 1*u 0*v ) 1e j 2π ( 1*u 1*v ) 1e j 2π (1*u 1*v ) 2e j 2π (1*u 0*v ) 1e j 2π (1*u 1*v ) j 2 sin 2πue j 2πv j 4 sin 2πu j 2 sin 2πue j 2πv j 4 sin 2πu (cos 2πv 1)Note: This signal is low pass in the horizontal direction v (weighted average),and band pass in the vertical direction u (difference). Yao Wang, 2016EL-GY 6123: Image and Video Processing37
Graph of F(u,v)2040608010020 40 60 80 100 Yao Wang, 2016du [-0.5:0.01:0.5];dv [-0.5:0.01:0.5];Fu abs(sin(2 * pi * du));Fv cos(2 * pi * dv);F 4 * Fu' * (Fv 1);mesh(du, dv, F);colorbar;Imagesc(F);colormap(gray); truesize;EL-GY 6123: Image and Video ProcessingUsing MATLAB freqz2:f [1,2,1;0,0,0;-1,-2,-1];freqz2(f)38
DSFT vs. 2D DFT 2D DSFT F(u, v) f (m, n)e j 2 π (mu nv)m n 2D DFTM 1 N 1F(k, l) f (m, n)e j 2 π (mk/M nl/N )m 0 n 0 (MxN) point 2D DFT are samples of DSFT for images ofsize MxN at u m/M, v n/N. Can be computed usingFFT algorithms. 1D FFT along rows, then 1D FFTalong columns Yao Wang, 2016EL-GY 6123: Image and Video Processing39
Display of DFT of ImagesImshow(img)MFimg abs(fft2(img)), imshow(MFimg,[ ])fftshift to shift the (0,0) freq.to centerLog mapping to enhancecontrastSMFimg fftshift(MFimg)imshow(SMFimg,[ ]) Yao Wang, 2016LSMFimg log(SMFimg 1);Imshow(LSMFimg,[ ])EL-GY 6123: Image and Video Processing40
DFT of Typical ImagesHow to interpret the DSFT image?Why is it bright at center and hassome line structures? Yao Wang, 2016 F(u,v) (obtained using 2D FFT)ff abs(fft2(f));imagesc(fftshift(log(ff 1)));Log mapping to enhance contrastFftshift to shift the (0,0) freq. to centerEL-GY 6123: Image and Video Processing41
Which one below is the DFT of which one above? Yao Wang, 2016EL-GY 6123: Image and Video Processing42
Yao Wang, 2016EL-GY 6123: Image and Video Processing43
DSFT of Separable Signals Separable signal:– h(m,n) hx(m) hy(n) 2D DSFT of separable signal product of 1D DSFT of each 1Dcomponent– H(u,v) Hx(u) Hy(v)– Hx(u): 1D FT of hx– Hy(v): 1D FT of hy Yao Wang, 2016EL-GY 6123: Image and Video Processing44
DSFT of Special Signals Constant - pulse at (0,0)Rectangle - 2D digital sincSinc - Rectangle (ideal low pass)Complex exponential with freq (u0,v0) - pulse at(u0,v0) Can be shown easily by making use of the fact that thesignal is separable Yao Wang, 2016EL-GY 6123: Image and Video Processing45
Using Separable Processing to ComputeDSFT 3x3 averaging filter! 1 1 1 ! 1 #&# &H 1 / 9 # 1 1 1 & 1 / 3# 1 &!" 1 1 1 %1 / 3 h1h1T# 1 1 1 &# 1 &"%" % Recognizing that the filter is separable1 / 3!" 1 1 1 # H1 (v) (1e j 2 π v 1 1e j 2 π v ) / 3 (1 2 cos2π v) / 3H (u, v) H1 (u)H1 (v) (1 cos2π u)(1 cos2π v) / 9 Yao Wang, 2016EL-GY 6123: Image and Video Processing46
Using Separable Processing to ComputeDSFT Another example 1 0 1 1 H 2 0 2 2 [1 0 1] hx hTy 1 0 1 1 Recognizing that the filter is separable[1[10 1] Fy (v) 1e j 2πv 0 ( 1)e j 2πv 2 j sin 2πv2 1] Fx (u ) 1e j 2πu 2 e j 2πu 2 2 cos 2πuF (u, v) Fx (u ) Fy (v) 4 j (1 cos 2πu ) sin 2πv Yao Wang, 2016EL-GY 6123: Image and Video Processing47
What about rotation? No theoretical proof, however, roughly it is still trueRotation in space - Rotation in freq.Example:f(m,n) δ(m) (horizontal line) - F(u,v) δ (v) (vertical line)What about f(m,n) δ(m-n) (diagonal line)?– - F(u,v) δ(u v) (antidiagonal line!)– Homework Yao Wang, 2016EL-GY 6123: Image and Video Processing48
Linear Convolution over Continuous Space 1D convolution( x ) h( x ) f Equalities f ( x α )h(α )dα f (α )h( x α )dαf ( x) δ ( x) f ( x),f ( x) δ ( x x0 ) f ( x x0 ) 2D convolutionf ( x , y ) h ( x, y ) Yao Wang, 2016f ( x α , y β )h(α , β )dαdβf (α , β )h(x α , y β )dαdβEL-GY 6123: Image and Video Processing49
Examples of 1D Convolutionf(x)f(α)h(x-α)h(x-α)10x1x-10αx(1) 0 x 1h(x)f(α)h(x-α)h(x-α)10αxx1x-1xα0αx-1 1(2) 1 x 2h(-α)f(x)*h(x)1101αx1 Yao Wang, 2016EL-GY 6123: Image and Video Processing2x50
Example of 2D Convolutiony1βf(x,y)1h(x-α,y-β)1yxx(1) 0 x 1, 0 y 1g(x,y) x*yy1h(x,y)1y2α1h(x-α,y-β)yx-111α(2) 0 x 1, 1 y 2g(x,y) x*(2-y)β1f(x,y)*h(x,y)αx(3) 1 x 2, 0 y 1g(x,y) (2-x)*y Yao Wang, 2016xβx2h(x-α,y-β) βy 1y-1yh(x-α,y-β)1y-11x-1 xα(4) 1 x 2, 1 y 2g(x,y) (2-x)*(2-y)xEL-GY 6123: Image and Video Processing51
Convolution of 1D Discrete Signals Definition of convolutionf ( n) * h( n) f ( n m) h ( m) m f ( m) h ( n m)m The convolution with h(n) can be considered as the weightedaverage in the neighborhood of f(n), with the filter coefficientsbeing the weights– sample f(n-m) is multiplied by h(m) The filter h(n) is the impulse response of the system (i.e. the outputto an input that is an impulse) Signal length before and after filtering– Original signal length: N– Filter length: K– Filtered signal length: N K-1 Yao Wang, 2016EL-GY 6123: Image and Video Processing52
Example of 1D Discrete Convolutionf(n)03h(n-m)nn-5h(n)05nn0(a) n 0, g(n) 0h(n-m)h(-m)0 nn-5-50mh(n-m)n Yao Wang, 2016m(b) 0 n 8, g(n) 0f(n)*h(n)0m8nn-50m(c) n 8, g(n) 0EL-GY 6123: Image and Video Processing53
Convolution of 2D Discrete Signals g(m, n) f (m, n)* h(m, n) f (m k, n l)h(k, l)k l f (k, l)h(m k, n l)k l Each new pixel g(m,n) is a weighted average of its neighboringpixels in the original image:– Pixel f(m-k,n-l) is weighted by h(k,l) We may use matrices to represent both signal (F) and filter (H) anduse F*H to denote the convolution Yao Wang, 2016EL-GY 6123: Image and Video Processing54
Example of 2D Discrete Convolutionmkf(k,l)h(-1-k, h(m,n)kk1h(-k,-l)12l Yao Wang, 2016lEL-GY 6123: Image and Video Processing355
Example: Averaging and Weighted Averaging1 9111111111100100 100100100100200 205203100100195 200200100100200 205195100100100 1001001001 16121242121100100 100100100100100 100100 100100144 167145100100156 176158 100100167 200168100100174 201175 100100144 166144100100156 175156 100100100 100100100100100 100100 100 Yao Wang, 2016EL-GY 6123: Image and Video Processing56
ExampleOriginal image Yao Wang, 2016Average filtered imageEL-GY 6123: Image and Video ProcessingWeighted Averagefiltered image57
What does h(m,n) mean? Any operation that is linear and shift invariant can bedescribed by a convolution with a filter h(m,n)! h(m,n) is the impulse response of the system (i.e.output of the system to an impulse input) Better known as point spread function, indicating how asingle point (i.e. an impulse) in the original image wouldbe spread out in the output image Yao Wang, 2016EL-GY 6123: Image and Video Processing58
Point Spread Function The point spread function of an imaging system (e.g. a camera or amedical imaging system) describes the resolution of the system:– Two object points cannot be separated if they are closer than the support of thepoint spread function!PSFXInput image Yao Wang, 2016Output imageEL-GY 6123: Image and Video Processing59
Boundary of Filtered Image An image of size MxN convolvingwith a filter of size KxL will yield animage of size (M K-1,N L-1) If the filter is symmetric with (2k 1)x(2k 1) samples, the convolvedimage should have an extraboundary of thickness k on eachside outside the original image(outer boundary). The values alongthe outer boundary depend on thepixel values outside the originalimage Filtered values in the innerboundary of k pixels inside theoriginal image also depend on thepixel values outside the originalimage Yao Wang, 2016Orange Red: original image sizeRed: Valid part of the output image(does not depend on pixels outsidethe original image)Orange: inner boundaryGreen: outer boundaryEL-GY 6123: Image and Video Processing60
Boundary Treatment: Zero Padding MxN image convolved with KxL filter - (M K-1)x(N L-1) imageFiltered values at the boundary of the output image depend on theassumed value of the pixels outside the original imageZero padding: Assuming pixel values are 0 outside the original image0000200 2050203000195 20020000200 20519500000Actual image pixels Yao Wang, 20160Extended pixels1 9111111111Outer boundary200/9(200 205)/9. (200 195)/9(200 205 195 200)/9 --- Inner boundaryEL-GY 6123: Image and Video Processing61
Boundary Treatment: Symmetric Extension Assuming pixels values outside the image are the same as their mirroringpixels inside the image– Lead to less discontinuity in the filtered image along the outer and innerboundaries200200200 205200 205203203203203195195 200200200200200 205195195200200 205195195Actual image pixels Yao Wang, 2016Extended pixels1 9111111111(200 200 200 200)/9(200 200 205 200 200 205)/9. (200 200 200 200 195 195)/9(200 200 205 200 200 205 195 195 200)/9 --- EL-GY 6123: Image and Video Processing62
Simplified Boundary Treatment Filtered image size original image sizeOnly compute in the valid region.– Assign 0 or keep original value in the inner boundary1 9200 205203111111111Outer boundary: not considered Yao Wang, 2016195 200200200 205195200205203195 200200205195Inner boundary: not processedEL-GY 6123: Image and Video Processing63
Example: Simplified Boundary Treatment1 9111111111100100 100100100100200 205203100100195 200200100100200 205195100100100 1001001001 16121242121100100 100100100100100 100100 100100144 167145100100156 176158 100100167 200168100100174 201175 100100144 166144100100156 175156 100100100 100100100100100 100100 100 Yao Wang, 2016EL-GY 6123: Image and Video Processing64
Sample Matlab Program (With SimplifiedBoundary Treatment)% readin bmp filex imread('lena.bmp');[xh xw] size(x);y double(x);% define 2D filterh ones(5,5)/25;[hh hw] size(h);hhh (hh - 1) / 2;hhw (hw- 1) / 2;% linear convolution, assuming the filter is non-separable (although this example filter is separable)z y; %or z zeros(xh,xw) if not low-pass filterfor m hhh 1:xh - hhh,%skip first and last hhh rows to avoid boundary problemsfor n hhw 1:xw - hhw,%skip first and last hhw columns to avoid boundary problemstmpy 0;for k -hhh:hhh,for l -hhw:hhw,tmpv tmpv y(m - k,n – l)* h(k hhh 1, l hhw 1);%h(0,0) is stored in h(hhh 1,hhw 1)endendz(m, n) tmpv;%for more efficient matlab coding, you can replace the above loop withz(m,n) sum(sum(y(m-hhh:m hhh,n-hhw:n hhw).*h))endend Yao Wang, 2016EL-GY 6123: Image and Video Processing65
Separable Filters A filter is separable if h(m, n) hx(m)hy(n). Matrix representationTH hx hy– Where hx and hy are column vectors Example 1 2 1 1 1 0 1 1 1111H x 0 0 0 0 [1 2 1]; H y 2 0 2 2 [ 1 0 1]4444 1 2 1 1 1 0 1 1 Yao Wang, 2016EL-GY 6123: Image and Video Processing66
Separable Filtering If h(m,n) is separable, the 2D convolution can be accomplished byfirst applying 1D filtering along each row using hy(n), and thenapplying 1D filtering to the intermediate result along each columnusing the filter hx(n) (or column filtering followed by row filtering) Prooff (m, n) * h(m, n) f (m k , n l )hx (k )hy (l )kl f (m k , n l )hy (l ) hx (k )k l g y (m k , n)hx (k )k ( f (m, n) * hy (n)) * hx (m) Yao Wang, 2016EL-GY 6123: Image and Video Processing67
Results Using Sobel FiltersOriginal imageFiltered image by HxFiltered image by Hy 1 2 1 1 1 0 1 1 1111H x 0 0 0 0 [1 2 1]; H y 2 0 2 2 [ 1 0 1]4444 1 2 1 1 1 0 1 1 What do Hx and Hy do? Yao Wang, 2016EL-GY 6123: Image and Video Processing68
Computation Cost: Non-Separable vs.Separable Filtering Suppose: Image size M*N, filter size K*L. Ignoring outer boundary pixelsNon-separable filtering– Weighted average on each pixel: K*L mul; K*L – 1 add.– For all pixels: M*N*K*L mul; M*N*(K*L-1) add.– When M N, K L: M2K2 mul M2(K2-1) add. Separable filtering:––––––––Each pixel in a row: L mul; L-1 add.Each row: N*L mul; N*(L-1) add.M rows: M*N*L mul; M*N*(L-1) add.Each pixel in a column: K mul; K-1 add.Each column: M*K mul; M*(K-1) add.N columns: N*M*K mul; N*M*(K-1) add.Total: M*N*(K L) mul; M*N*(K L-2) add.When M N, K L: 2M2K mul; 2M2(K-1) add. Significant savings if K (and L) is large! Yao Wang, 2016EL-GY 6123: Image and Video Processing69
MATLAB Function: conv2( ) C conv2(H,F,shape)– Shape ‘full’ (Default): C includes both outer and innerboundary, using zero padding– Shape “same”: C includes the inner boundary, using zeropadding– Shape “valid’: C includes the convolved image without theinner boundary, computed without using pixels outside theoriginal image C conv2(h1,h2,F,shape)– Separable filtering with h1 for column filtering and h2 for rowfiltering Yao Wang, 2016EL-GY 6123: Image and Video Processing70
Notes about MATLAB implementation Input image needs to be converted to integer or double– Never do numerical operations on unsigned character type! Output image value may not be in the range of (0,255) and may not be integersTo display or save the output image properly– Renormalize to (0,255) using a two-pass operation First pass: save directly filtered value in an intermediate floating-point array Second pass: find minimum and maximum values of the intermediateimage, renormalize to (0,255) and rounding to integers– F round((F1-fmin)*255/(fmax-fmin))– To display the unnormalized image directly in MATLAB, useimagesc(img). Or imshow(img, [ ]). Yao Wang, 2016EL-GY 6123: Image and Video Processing71
Convolution Theorem Convolution Theoremf *h F H, Prooff h F *Hg (m, n) f (m, n) * h(m, n) f (m k , n l )h(k , l )klFT on both sidesG (u , v) f (m k , n l )h(k , l )e j 2π ( mu nv )m , n k ,l f (m k , n l )e j 2π (( m k ) u ( n l ) v ) h(k , l )e j 2π ( ku lv )m , n k ,l f (m k , n l )e j 2π (( m k ) u ( n l ) v ) h(k , l )e j 2π ( ku lv )m,nk ,l f (m' , n' )e j 2π ( m 'u n 'v ) h(k , l )e j 2π ( ku lv )m ', n 'k ,l F (u , v) H (u , v) Yao Wang, 2016EL-GY 6123: Image and Video Processing72
Another view of convolution theoremf (m, n)* h(m, n) F(u, v)H (u, v) F(u,v)H(u,v) Modifying the signal’s each frequency component’complex magnitude F(u,v) by H(u,v)1/ 2f (m, n) 1/ 2 1/ 2 1/ 2g(m, n) F (u, v)e1/21/2 1/2 1/2 j 2π ( mu nv )dudvF(u, v)H (u, v)e j 2 π (mu nv) du dvH(u,v) is also called Frequency Response of the 2D LSI system––– Yao Wang, 2016exp{2π (um vn)} - H(u,v) exp{2π (um vn)} H(u,v) exp{2π (um vn) P(H(u,v))}Sinusoid (or complex exponential) input - sinusoid (complex exponential) output!H(u,v) describes how the magnitude and phase of a sinusoid input with frequency (u,v) arechanged!EL-GY 6123: Image and Video Processing73
Explanation of Convolutionin the Frequency Domainh(x)f(x)-Δ/2Δ/2x-Δ/4F(u)-1/Δ Yao Wang, 2016g(x) f(x)*h(x)Δ/4xxG(u) F(u)H(u)H(u)1/Δu-2/Δ2/ΔuEL-GY 6123: Image and Video Processing-1/Δ1/Δu74
Example Given a 2D filter, determine its frequency response. Apply to agiven image, show original image and filtered image in pixel andfreq. domain 1 11 1h 25 1 1 Yao Wang, 20161 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 EL-GY 6123: Image and Video Processing75
1 11 1h 25 1 1 Yao Wang, 20161 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 EL-GY 6123: Image and Video Processing76
Matlab Program Usedx imread('lena256.bmp');figure(1); imshow(x);f double(x);ff abs(fft2(f));figure(2); imagesc(fftshift(log(ff 1))); colormap(gray);truesize;axis off;h ones(5,5)/9;hf abs(freqz2(h));figure(3);imagesc((log(hf 1)));colormap(gray);truesize;axis off;y conv2(f, is off;yf abs(fft2(y));figure(5);imagesc(fftshift(log(yf 1)));colormap(gray);truesize;axis off; Yao Wang, 2016EL-GY 6123: Image and Video Processing77
1 1 1 1H1 0 0 0 3 1 1 1 Yao Wang, 2016EL-GY 6123: Image and Video Processing78
Typical Filter TypesLow PassHigh Passuu0.5-0.5u0.50.5-0.5Band Passv-0.50.50.5v-0.5-0.50.5v-0.5Non-zero frequency components, where F(u,v) ! 0 Yao Wang, 2016EL-GY 6123: Image and Video Processing79
Calculate Linear Convolution Using DFT 1D case– f(n) is length N1, h(n) is length N2– g(n) f(n)*h(h) is length N N1 N2-1.– To use DFT, need to extend f(n) and h(n) to length N by zeropadding.– H(k) can be precalculatedf(n)(zero padded)*N-point DFTDFTF(k) Yao Wang, licationEL-GY 6123: Image and Video ProcessingG(k)80
Computing 2D Convolution Using 2D DFTZeropaddingRelation between spatial and frequency domain operation:g ( x, y ) h( x, y ) f ( x, y ) G (u, v) H (u, v) F (u, v)h( x, y ) IDFT ( H (u, v)), H (u, v) DFT (h( x, y )).Typically DFT size image size. This corresponds to circular convolution, which differsfrom linear convolution at the inner boundaries. Only correct in the valid region. Yao Wang, 2016EL-GY 6123: Image and Video Processing81
Image Filtering Using DFT Typically DFT size
Two-Dimensional Fourier Transform and Linear Filtering Yao Wang . Image and Video Processing 14 Two Dimension Continuous Space Fourier Transform (CSFT) Basis functions Forward – Transform . – For separable signal, one can simply compute two 1D transforms and take their product! F 2 {f (x, y)} F y {F x
The Matlab Hilbert transform operates on one-dimensional data. For the work described here, it is necessary to adapt the HT to two-dimensional data. I do this by analogy to the two-dimensional Fourier Transform (strictly, the "discrete Fourier Transform"). Basic equations (e.g. Lim, 1990) show that the two-dimensional discrete Fourier
to denote the Fourier transform of ! with respect to its first variable, the Fourier transform of ! with respect to its second variable, and the two-dimensional Fourier transform of !. Variables in the spatial domain are represented by small letters and in the Fourier domain by capital letters. expressions, k is an index assuming the two values O
(d) Fourier transform in the complex domain (for those who took "Complex Variables") is discussed in Appendix 5.2.5. (e) Fourier Series interpreted as Discrete Fourier transform are discussed in Appendix 5.2.5. 5.1.3 cos- and sin-Fourier transform and integral Applying the same arguments as in Section 4.5 we can rewrite formulae (5.1.8 .
FT Fourier Transform DFT Discrete Fourier Transform FFT Fast Fourier Transform WT Wavelet Transform . CDDWT Complex Double Density Wavelet Transform PCWT Projection based Complex Wavelet Transform viii. . Appendix B 150 Appendix C 152 References 153 xiii.
straightforward. The Fourier transform and inverse Fourier transform formulas for functions f: Rn!C are given by f ( ) Z Rn f(x)e ix dx; 2Rn; f(x) (2ˇ) n Z Rn f ( )eix d ; x2Rn: Like in the case of Fourier series, also the Fourier transform can be de ned on a large class of generalized functions (the space of tempered
Expression (1.2.2) is called the Fourier integral or Fourier transform of f. Expression (1.2.1) is called the inverse Fourier integral for f. The Plancherel identity suggests that the Fourier transform is a one-to-one norm preserving map of the Hilbert space L2[1 ;1] onto itself (or to anoth
option price and for the Fourier transform of the time value of an option. Both Fourier transforms are expressed in terms of the characteristic function of the log price. 3.1 . The Fourier Transform of an Option Price Let k denote the log of the strike price K, and let C T (k) be the desired value of a T-maturity call option with strike exp(k
1 eng1a01 1 transactions essential english language skills 4 3 7 2 eng1a02 1 ways with words literatures in english 5 3 9 3 eng2a03 2 writing for academic and professional 4 4 11 . 3 success 4 eng2a04 2 zeitgeist readings on contempo rary culture 5 4 13 5 eng3a05 3 signatures expressing the self 5 4 15 6 eng4a06 4 spectrum literature and contemporary issues 5 4 17 to tal 22 .
