Matlab Script implementing MACH Filter - Matlab code for Processing Image | MAP 6938, Study Guides, Projects, Research of Mathematics

Download Matlab Script implementing MACH Filter - Matlab code for Processing Image | MAP 6938 and more Study Guides, Projects, Research Mathematics in PDF only on Docsity! Subhabrata Bhattacharya/Brian Brennan MAP 6938 PROJECT 02 0.1 Theory 1. Use Newton’s method to give a numerical approach for solving: ∇f(x, y) = 0 The Newton’s method for solving g(x) = 0, is given by : g(x) ≈ g(a) + g′(a)(x− a) =0 [where a is an initial guess] or, x =a− g(a) g′(a) This can be interpolated to the following for nth value of x : or,xn+1 = xn − g(xn) g′(xn) [For x0 = a] Since our function comprises of two variables x, y let us suppose : ~r = [ x y ] and let ~r0 = [ a b ] be our initial guess. Then : ∇f(x, y) ≈ ∇f(a, b) +∇(∇f(a, b))(~r − ~r0) = 0 or,∇f(x, y) ≈ ∇f(a, b) +∇2f(a, b)(~r − ~r0) = 0 Now, ∇2f(x, y) = [ ∂2f ∂x2 ∂2f ∂x∂y ∂2f ∂y∂x ∂2f ∂y2 ] = Hf(x, y) which is the Hessian of f(x, y). From the above we could conclude that : ∇f(x, y) ≈ ∇f(a, b) +Hf(x, y)(~r − ~r0) =0 or, ~r =~r0 −Hf(a, b)−1∇f(a, b) (1) We could generalize from (1) for Xn and Yn :[ xn+1 yn+1 ] = [ xn yn ] − [Hf(xn, yn)]−1∇f(xn, yn) [For [ x0 y0 ] = [ a b ] ] (2) 2. Compare the formulas from the previous problem with those for the steepest descent methods of finding the minimum of f(x, y). Derive any formula that is needed. The minimum of f(x, y) can be approximated by iterating the following sequence:[ Xn+1 Y n+ 1 ] = [ Xn Y n ] − αn∇f(Xn, Yn) where [ X0 Y0 ] = [ a b ] Comparing this sequence to (2), it is clear that they are equivalent when αn = [Hf(Xn, Yn)] −1. In the previous problem we were using Newton’s method to find a zero of the gradient. This of course will lead us to a stationary point which would be a minimum value. So it should be expected that Newton’s method, applied to the gradient would be related to a method for approximating a minimum. This method is the Steepest Descent. 3. Let A be a nxn symmetric matrix and b be a column vector of dimension n. Derive the following: (a) ∇f(~x) if f(~x) = bT~x For the sake of simplicity, let us assume: ~b = [ b1 b2 ] , ~x = [ x1 x2 ] (3) Therefore, f(~x) = ~bT~x = b1x1 + b2x2 (4) 1 Subhabrata Bhattacharya/Brian Brennan MAP 6938 PROJECT 02 Now, ∇f(~x) = [ ∂f(~x) ∂x1 ∂f(~x) ∂x2 ] Taking partial derivatives on (4) wrt x1, x2 we get: ∇f(~x) = [ b1 b2 ] = ~bT We could easily infer from the above that for ~b =  b1 b2 ... bn  and ~x =  x1 x2 ... xn  The gradient ∇f(~x) is given by :  ∂f(~x) ∂x1 ∂f(~x) ∂x2 ... ∂f(~x) ∂xn  = [ b1 b2 . . . bn ] = ~bT (5) (b) ∇f(~x) if f(~x) = ~xTA~x Let us suppose A = [ a1,1 a1,2 a2,1 a2,2 ] Therefore using (3) we have, f(~x) = ~xTA~x = [ x1 x2 ] [ a1,1 a1,2 a2,1 a2,2 ] [ x1 x2 ] = a1,1x12 + (a1,2 + a2,1)x1x2 + a2,2x22 (6) The gradient ∇f(~x) is given by :[ ∂f(~x) ∂x1 ∂f(~x) ∂x2 ] = [ 2a1,1x1 (a1,2 + a2,1)x2 (a1,2 + a2,1)x1 2a2,2x2 ] = [ 2a1,1 (a1,2 + a2,1) (a1,2 + a2,1) 2a2,2 ] [ x1 x2 ] = [ A+AT ] [ x1 x2 ] = 2A~x [As A is symmetric, AT = A] (7) The above relation could be generalized an nxn symmetric matrix A and n dimensional column vectors ~x and ~b. (c) Find the Hessian of g(x) = XTAx− bTx The Hessian of g(x), H(g)i,j(~x) = ∂ ∂xi ∂ ∂xj g(~x) = ∇2g(~x) However from (5) and (7), we could conclude : ∇g(~x) = 2A~x− bT So, ∇2g(~x) = 2A Thus, Hg(~x) = 2A 4. Let f and g be functions of two variables. Derive a formula for the discrete Fourier transform of f ⊗ g in terms of those for f and g. Let us assume, f(x, y) and g(x + p, y + q) be our bi-variate functions. Therefore the correlation function between these functions is given by : f(x, y)⊗ g(x+ p, v + q) = ∑ p ∑ q f(x, y)g(x+ p, y + q) (8) 2 Subhabrata Bhattacharya/Brian Brennan MAP 6938 PROJECT 02 42 43 % Compute the Average S im i l a r i t y measure exp r e s s i on 44 SPX = ze ro s ( nPix , 1 ) ; 45 f o r l c n t = 1 : nImg 46 PXdif f = PX( : , l c n t ) − PXavg ; 47 SPX = SPX + conj ( PXdif f ) . *PXdif f ; 48 end 49 SPX = SPX/nImg ; 50 % Compute h f o r each s e t o f t r a i n i n g images . h i s hence a nPix x [ number o f 51 % s e t s o f t r a i n i n g images ] 52 h ( : , kcnt ) = (1 . /(SPX + ∆* ones ( nPix , 1 ) ) ) . * xavg ; 53 kcnt = kcnt +1; 54 j cn t = j cn t + 45 ; 55 end 56 end 57 %Number o f f i l t e r s 58 n f i l t = kcnt −1; 59 60 % Test ing i f the code works 61 % Given an image : f i nd out the nea r e s t one f i t e r c l a s s that has 62 % c o r r e l a t i o n s output s im i l a r to i t 63 64 fnames = ' /home/subh/ cour s e s /map6398/ p r o j e c t s / pro j e c t −2/datase t /* . i s 4 ' ; 65 data = ' /home/subh/ cour s e s /map6398/ p r o j e c t s / pro j e c t −2/datase t / ' ; 66 row = 128 ; 67 c o l = 128 ; 68 d = row* c o l ; 69 70 % Pick an image randomly from the datase t 71 f i l e s = d i r ( fnames ) ; 72 f l i d x = f i nd (¬ [ f i l e s . i s d i r ] ) ; 73 nImg = length ( f l i d x ) ; 74 idx = length ( f l i d x ) − round (nImg* rand ( 1 ) ) ; 75 76 flname = f i l e s ( f l i d x ( idx ) ) .name ; 77 d i sp ( f lname ) ; 78 79 % Load image 80 fd = fopen ( s t r c a t ( data , f lname ) , ' r ' , ' l ' ) ; 81 img = f r ead ( fd , [ 1 512*512 ] , ' ubit8 ' ) ; 82 f c l o s e ( fd ) ; 83 img = img (129 : end ) ; 84 img = reshape ( img , [100 1 0 0 ] ) ; 85 img = imre s i z e ( img , [ row co l ] , ' b i cub i c ' ) ; 86 Xtest = reshape ( img , d , 1 ) ; 87 Xtest = angle ( f f t 2 ( Xtest ) ) ; 88 % Computing the r e sponse s o f the amplitude normal ized t e s t image in 89 % frequency domain with the s e t s o f f i l t e r 90 91 f o r i c n t =1: n f i l t 92 g ( : , i c n t ) = Xtest . *h ( : , i c n t ) ; 93 resp ( : , i c n t ) = sum( g ( : , i c n t ) ) ; 94 end 95 Cidx = f i nd ( resp==max( resp ) ) ; 96 % In t e r p r e t i n g the i n d i c e s 97 % This i s the o rgan i z a t i on o f c l a s s e s : f o r each e l e v a t i o n the re are 8 s e t s 98 % of f i l t e r Cidx/8 would g ive us the e l e v a t i o n c l a s s the f i l t e r b e l ong s . 99 eClas s = 1 + ( Cidx /8)*3 ; 100 aClass = 45*mod(Cidx , 8 ) ; 101 d i sp ( eClass ) ; 102 d i sp ( aClass ) ; Listing 1: Matlab script implementing MACH filter 5 Subhabrata Bhattacharya/Brian Brennan MAP 6938 PROJECT 02 1 func t i on [ spatM , freqM ] = get InputMatr i ce s ( eRange , aRange ) 2 % This func t i on loads f i l e s with g iven input e l e v a t i o n and output Ranges. 3 % Returns spatM , freqM matr i ce s which are s p a t i a l domain and frequency 4 % domain trans forms o f the images in the s e l e c t e d e l e v a t i o n and azimuth 5 % ranges 6 d i rdata = ' /home/subh/ cour s e s /map6398/ p r o j e c t s / pro j e c t −2/datase t / ' ; 7 row = 128 ; 8 c o l = 128 ; 9 d = row* c o l ; 10 11 % Al l images in the datase t are o f the f i l e t y p e Exxaxxx. i s4 where x i s an 12 % i n t e g e r . Also xx and xxx are e i t h e r 0 or mu l t i p l e s o f t h r e e . 13 kcnt = 1 ; 14 f o r i c n t = eRange ( 1 ) : 3 : eRange (2 ) 15 f o r j cn t = aRange ( 1 ) : 3 : aRange (2 ) 16 i f ( i c n t < 10) 17 e l e v = [ ' 0 ' num2str ( i c n t ) ] ; 18 e l s e 19 e l e v = num2str ( i c n t ) ; 20 end 21 22 i f i c n t == 0 23 i f ( j cn t < 10) 24 a z i = [ ' 00 ' num2str ( j cn t ) ] ; 25 e l s e i f ( j cn t < 100) 26 a z i = [ ' 0 ' num2str ( j cn t ) ] ; 27 e l s e 28 a z i = num2str ( j cn t ) ; 29 end 30 e l s e 31 a z i = num2str ( j cn t ) ; 32 end 33 34 fname = [ d i rdata 'E ' e l e v ' a ' a z i ' . i s 4 ' ] ; 35 fd = fopen ( fname , ' r ' ) ; 36 i f fd == −1 37 cont inue ; 38 end 39 img = f r ead ( fd , [ 1 512*512 ] , ' ubit8 ' ) ; 40 f c l o s e ( fd ) ; 41 img = img (129 : end ) ; 42 img = reshape ( img , [100 1 0 0 ] ) ' ; 43 img = imre s i z e ( img , [ row co l ] , ' b i cub i c ' ) ; 44 45 % Test ing i f the image i s loaded 46 % f i g u r e ( 1 ) , colormap ( gray (256 ) ) , imagesc ( img ) , ax i s ( ' image ' ) , drawnow 47 % Arranging a l l the images as a c o l l e c t i o n o f column vec to r s spatM 48 spatM ( : , kcnt ) = reshape ( img , d , 1 ) ; 49 % Arranging a l l the Frequency domain trans form o f the images as a 50 % c o l l e c t i o n o f column vec to r s freqM 51 freqM ( : , kcnt ) = reshape ( f f t 2 ( img ) , d , 1 ) ; 52 kcnt = kcnt + 1 ; 53 end 54 end Listing 2: Matlab code for processing Image 6 Subhabrata Bhattacharya/Brian Brennan MAP 6938 PROJECT 02 0.2.2 Results Training image Filename Identified Elevation Identified Azimuth E06a258.is4 6.2500 270 E30a279.is4 53.5000 180 E84a135.is4 52 0 E63a336.is4 64.1250 335 E39a105.is4 41.2500 70 E03a168.is4 44.1250 135 E03a0.is4 8 0 E78a225.is4 80.1250 20 E00a300.is4 04.1250 235 E03a174.is4 5 180 0.2.3 Areas of Improvements The results are not entirely as expected. The deviations from the expected results could mainly be due to the unavailability of training images for certain images belonging to a particular set of orientations. 7