Digital Image Processing_S. Jayaraman, S. Esakkirajan and T. Veerakumar

Scilab Textbook Companion for Digital Image Processing by S. Jayaraman, S. Esakkirajan And T. Veerakumar1 Created by R.Senthilkumar M.E., Assistant Professor Electronics Engineering Inst. of Road and Transport Tech., Erode College Teacher Dr. R.K.Gnanamurhty, Principal,V. C. E. W. Cross-Checked by Anuradha Amrutkar, IIT Bombay August 10, 2013

1 Funded

by a grant from the National Mission on Education through ICT, http://spoken-tutorial.org/NMEICT-Intro. This Textbook Companion and Scilab codes written in it can be downloaded from the ”Textbook Companion Project” section at the website http://scilab.in

Book Description Title: Digital Image Processing Author: S. Jayaraman, S. Esakkirajan And T. Veerakumar Publisher: Tata McGraw - Hill Education Pvt. Ltd, New Delhi Edition: 3 Year: 2010 ISBN: 978-0-07-014479-8

1

Scilab numbering policy used in this document and the relation to the above book. Exa Example (Solved example) Eqn Equation (Particular equation of the above book) AP Appendix to Example(Scilab Code that is an Appednix to a particular Example of the above book) For example, Exa 3.51 means solved example 3.51 of this book. Sec 2.3 means a scilab code whose theory is explained in Section 2.3 of the book.

2

Contents List of Scilab Codes

4

1 Introduction to Image Processing System

9

2 2D Signals and Systems

12

3 Convolution and Correlation

15

4 Image Transforms

24

5 Image Enhancement

34

6 Image Restoration and Denoising

44

7 Image Segmentation

64

8 Object Recognition

73

9 Image Compression

77

10 Binary Image Processing

81

11 Colur Image Processing

85

12 Wavelet based Image Processing

99

3

List of Scilab Codes Exa 1.3 Exa Exa Exa Exa Exa Exa Exa Exa Exa

1.13 2.12 2.16 3.1 3.2 3.3 3.7 3.8 3.11

Exa 3.12 Exa 3.13 Exa Exa Exa Exa Exa Exa Exa Exa

3.14 3.15 3.16 3.17 3.18 4.4 4.5 4.6

Exa 4.10 Exa 4.12 Exa 4.13

Program to calculate number of samples required for an image . . . . . . . . . . . . . . . . . . . . . . . . . . . False contouring Scilab code . . . . . . . . . . . . . . . Frequency Response . . . . . . . . . . . . . . . . . . . Frequency Response . . . . . . . . . . . . . . . . . . . 2D Linear Convlolution . . . . . . . . . . . . . . . . . 2D Linear Convolution . . . . . . . . . . . . . . . . . . 2D Linear Convolution . . . . . . . . . . . . . . . . . . 2D Linear Convolution . . . . . . . . . . . . . . . . . . 2D Linear Convolution . . . . . . . . . . . . . . . . . . Linear Convolution of any signal with an impule signal given rise to the same signal . . . . . . . . . . . . . . . Circular Convolution between two 2D matrices . . . . Circular Convolution exspressed as linear convolution plus alias . . . . . . . . . . . . . . . . . . . . . . . . . Linear Cross correlation of a 2D matrix . . . . . . . . Circular correlation between two signals . . . . . . . . Circular correlation between two signals . . . . . . . . Linear auto correlation of a 2D matrix . . . . . . . . . Linear Cross correlation of a 2D matrix . . . . . . . . DFT of 4x4 grayscale image . . . . . . . . . . . . . . . 2D DFT of 4X4 grayscale image . . . . . . . . . . . . Scilab code to intergchange phase information between two images . . . . . . . . . . . . . . . . . . . . . . . . Program to compute discrete cosine transform . . . . . Program to perform KL tranform for the given 2D matrix Program to find the singular value decomposition of given matrix . . . . . . . . . . . . . . . . . . . . . . . 4

9 10 12 12 15 15 16 17 17 18 18 19 20 20 21 22 22 24 24 25 26 28 32

Exa Exa Exa Exa Exa Exa Exa Exa Exa Exa Exa

5.5 5.7 5.9 5.13 5.16 5.20 6.1 6.5 6.7 6.9 6.13

Exa Exa Exa Exa Exa Exa Exa

6.18 6.21 6.23 6.24 7.23 7.25 7.27

Exa Exa Exa Exa

7.30 8.4 8.5 8.9

Exa Exa Exa Exa

9.9 9.59 10.17 10.19

Exa 11.4 Exa 11.12 Exa 11.16 Exa 11.18 Exa 11.21

Scilab code for brightness enhancement . . . . . . . . Scilab code for brightness suppression . . . . . . . . . Scilab code for Contrast Manipulation . . . . . . . . . Scilab code to determine image negative . . . . . . . . Scilab code that performs threshold operation . . . . . Program performs gray level slicing without background Scilab code to create motion blur . . . . . . . . . . . . Scilab code performs inverse filtering . . . . . . . . . . Scilab code performs inverse filtering . . . . . . . . . . Scilab code performs Pseudo inverse filtering . . . . . Scilab code to perform wiener filtering of the corrupted image . . . . . . . . . . . . . . . . . . . . . . . . . . . Scilab code to Perform Average Filtering operation . . Scilab code to Perform median filtering . . . . . . . . Scilab code to Perform median filtering of colour image Scilab code to Perform Trimmed Average Filter . . . . Scilab code for Differentiation of Gaussian function . . Scilab code for Differentiation of Gaussian Filter function Scilab code for Edge Detection using Different Edge detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . Scilab code to perform watershed transform . . . . . . To verify the given matrix is a covaraince matrix . . . To compute the covariance of the given 2D data . . . Develop a perceptron AND function with bipolar inputs and targets . . . . . . . . . . . . . . . . . . . . . . . . Program performs Block Truncation Coding BTC . . . Program performs Block Truncation Coding . . . . . . Scilab Code for dilation and erosion process . . . . . . Scilab Code to perform an opening and closing operation on the image . . . . . . . . . . . . . . . . . . . . Read an RGB image and extract the three colour components red green blue . . . . . . . . . . . . . . . . . . Read a Colour image and separate the colour image into red green and blue planes . . . . . . . . . . . . . . . . Compute the histogram of the colour image . . . . . . Perform histogram equalisation of the given RGB image This program performs median filtering of the colour image . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

34 36 36 38 38 40 44 45 48 50 54 55 57 59 61 64 67 69 71 73 73 75 77 79 81 82 85 86 88 90 91

Exa Exa Exa Exa

11.24 11.28 11.30 11.32

Fitlering only the luminance component . . . . . . . . Perform gamma correction for the given colour image . Perform Pseudo Colouring Operation . . . . . . . . . . Read an RGB image and segment it using the threshold method . . . . . . . . . . . . . . . . . . . . . . . . . . Exa 12.9 Scilab code to perform wavelet decomposition . . . . . Exa 12.42 Scilab code to generate different levels of a Gaussian pyramid . . . . . . . . . . . . . . . . . . . . . . . . . . Exa 12.57 Scilab code to implement watermarking in spatial domain Exa 12.63 Scilab code to implement wavelet based watermarking AP 1 2D Fast Fourier Transform . . . . . . . . . . . . . . . AP 2 2D Inverse FFT . . . . . . . . . . . . . . . . . . . . . AP 3 Median Filtering function . . . . . . . . . . . . . . . . AP 4 To caculate gray level . . . . . . . . . . . . . . . . . . AP 5 To change the gray level of gray image . . . . . . . . .

6

92 93 95 96 99 100 101 104 107 108 108 109 109

List of Figures 1.1

False contouring Scilab code . . . . . . . . . . . . . . . . . .

11

2.1 2.2

Frequency Response . . . . . . . . . . . . . . . . . . . . . . Frequency Response . . . . . . . . . . . . . . . . . . . . . .

13 14

4.1 4.2 4.3

Scilab code to intergchange phase information between two images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Program to compute discrete cosine transform . . . . . . . . Program to compute discrete cosine transform . . . . . . . .

27 29 30

5.1 5.2 5.3 5.4 5.5 5.6

Scilab code for brightness enhancement . . . . . . . . . . Scilab code for brightness suppression . . . . . . . . . . . Scilab code for Contrast Manipulation . . . . . . . . . . Scilab code to determine image negative . . . . . . . . . Scilab code that performs threshold operation . . . . . . Program performs gray level slicing without background

35 37 39 39 41 43

6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10

Scilab Scilab Scilab Scilab Scilab Scilab Scilab Scilab Scilab Scilab

7.1 7.2

Scilab code for Differentiation of Gaussian function . . . . . Scilab code for Differentiation of Gaussian function . . . . .

code code code code code code code code code code

. . . . . .

. . . . . .

to create motion blur . . . . . . . . . . . . . . . performs inverse filtering . . . . . . . . . . . . . performs inverse filtering . . . . . . . . . . . . . performs inverse filtering . . . . . . . . . . . . . performs Pseudo inverse filtering . . . . . . . . . to perform wiener filtering of the corrupted image to Perform Average Filtering operation . . . . . to Perform median filtering . . . . . . . . . . . . to Perform median filtering of colour image . . . to Perform Trimmed Average Filter . . . . . . .

7

45 47 49 51 53 56 58 60 61 63 65 66

7.3 7.4 7.5 7.6

Scilab Scilab Scilab Scilab

code code code code

for Differentiation of Gaussian Filter function . . for Differentiation of Gaussian Filter function . . for Edge Detection using Different Edge detectors to perform watershed transform . . . . . . . . .

9.1

Program performs Block Truncation Coding . . . . . . . . .

80

10.1 Scilab Code for dilation and erosion process . . . . . . . . . 10.2 Scilab Code to perform an opening and closing operation on the image . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

11.1 Read an RGB image and extract the three colour components red green blue . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Read a Colour image and separate the colour image into red green and blue planes . . . . . . . . . . . . . . . . . . . . . . 11.3 This program performs median filtering of the colour image . 11.4 Fitlering only the luminance component . . . . . . . . . . . 11.5 Perform gamma correction for the given colour image . . . . 11.6 Perform Pseudo Colouring Operation . . . . . . . . . . . . . 11.7 Read an RGB image and segment it using the threshold method

67 68 70 72

84 87 89 92 94 95 97 98

12.1 Scilab code to generate different levels of a Gaussian pyramid 102 12.2 Scilab code to implement watermarking in spatial domain . . 104

8

Chapter 1 Introduction to Image Processing System

Scilab code Exa 1.3 Program to calculate number of samples required for an image 1 2 3 4 5 6 7 8 9 10 11 12 13 14

// C a p t i o n : Program t o c a l c u l a t e number o f s a m p l e s r e q u i r e d f o r an image // Example1 . 3 // p a g e 12 clc ; close ; // d i m e n s i o n o f t h e image i n i n c h e s m = 4; n = 6; N = 400; // number o f d o t s p e r i n c h i n e a c h d i r e c t i o n N2 = 2* N ; // number o f d o t s p e r i n c h i n b o t h horizontal & vertical Fs = m * N2 * n * N2 ; disp ( Fs , ’ Number o f s a m p l e s r e u q i r e d t o p r e s e r v e t h e i n f o r m a t i o n i n t h e image= ’ ) // R e s u l t // Number o f s a m p l e s r e u q i r e d t o p r e s e r v e t h e i n f o r m a t i o n i n t h e image= 9

15

// 1 5 3 6 0 0 0 0 . check Appendix AP 4 for dependency: gray.sci check Appendix AP 5 for dependency: grayslice.sci

Scilab code Exa 1.13 False contouring Scilab code 1 // C a p t i o n : F a l s e c o n t o u r i n g S c i l a b c o d e 2 // F i g 1 . 1 3 3 // p a g e 13 4 clc ; 5 close ; 6 a = ReadImage ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 \ t i g e r p u b . j p g

’ ); 7 a = uint8 ( a ) ; 8 figure 9 imshow ( a ) 10 title ( ’ O r i g i n a l image ’ ) ; 11 // u s i n g 128 g r a y l e v e l s 12 figure 13 a_128 = grayslice (a ,128) ; 14 gray_128 = gray (128) ; 15 ShowImage ( a_128 , ’ Image w i t h 128 g r a y 16 17 18 19 20 21 22 23

levels ’,

gray_128 ) ; // u s i n g 64 g r a y l e v e l s figure a_64 = grayslice (a ,64) ; gray_64 = gray (64) ; ShowImage ( a_64 , ’ Image w i t h 64 g r a y l e v e l s ’ , gray_64 ) ; // u s i n g 32 g r a y l e v e l s figure a_32 = grayslice (a ,32) ; 10

Figure 1.1: False contouring Scilab code 24 gray_32 = gray (32) ; 25 ShowImage ( a_32 , ’ Image w i t h 32 g r a y l e v e l s ’ , gray_32 ) ; 26 // u s i n g 16 g r a y l e v e l s 27 figure 28 a_16 = grayslice (a ,16) ; 29 gray_16 = gray (16) ; 30 ShowImage ( a_16 , ’ Image w i t h 16 g r a y l e v e l s ’ , gray_16 ) ; 31 // u s i n g 8 g r a y l e v e l s 32 a_8 = grayslice (a ,8) ; 33 gray_8 = gray (8) ; 34 ShowImage ( a_8 , ’ Image w i t h 8 g r a y l e v e l s ’ , gray_8 ) ;

11

Chapter 2 2D Signals and Systems

Scilab code Exa 2.12 Frequency Response 1 // C a p t i o n : F r e q u e n c y R e s p o n s e 2 // F i g 2 . 1 2 3 // p a g e 60 4 clc ; 5 close ; 6 [X , Y ] = meshgrid ( - %pi :.09: %pi ) ; 7 Z = 2* cos ( X ) +2* cos ( Y ) ; 8 surf (X ,Y , Z ) ; 9 xgrid (1)

Scilab code Exa 2.16 Frequency Response 1 2 3

// C a p t i o n : F r e q u e n c y R e s p o n s e // F i g 2 . 1 6 // p a g e 64 12

Figure 2.1: Frequency Response 4 clc ; 5 close ; 6 [X , Y ] = meshgrid ( - %pi :.05: %pi ) ; 7 Z = 2 - cos ( X ) - cos ( Y ) ; 8 surf (X ,Y , Z ) ; 9 xgrid (1)

13

Figure 2.2: Frequency Response

14

Chapter 3 Convolution and Correlation

Scilab code Exa 3.1 2D Linear Convlolution 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

// C a p t i o n : 2−D L i n e a r C o n v o l u t i o n // Example3 . 1 & Example3 . 4 // p a g e 85 & p a g e 107 clc ; x =[4 ,5 ,6;7 ,8 ,9]; h = [1;1;1]; disp (x , ’ x= ’ ) disp (h , ’ h= ’ ) [y ,X , H ] = conv2d2 (x , h ) ; disp (y , ’ L i n e a r 2D c o n v o l u t i o n r e s u l t y = ’ ) // R e s u l t // L i n e a r 2D c o n v o l u t i o n r e s u l t y = // // 4. 5. 6. // 11. 13. 15. // 11. 13. 15. // 7. 8. 9.

Scilab code Exa 3.2 2D Linear Convolution 15

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

// C a p t i o n : 2−D L i n e a r C o n v o l u t i o n // Example3 . 2 & Example3 . 5 & Example3 . 9 // p a g e 91 & p a g e 108 & p a g e 116 clc ; x =[1 ,2 ,3;4 ,5 ,6;7 ,8 ,9]; h = [1 ,1;1 ,1;1 ,1]; y = conv2d2 (x , h ) ; disp (y , ’ L i n e a r 2D c o n v o l u t i o n r e s u l t y = ’ ) // R e s u l t // L i n e a r 2D c o n v o l u t i o n r e s u l t y = // // 1. 3. 5. 3. // 5. 12. 16. 9. // 12. 27. 33. 18. // 11. 24. 28. 15. // 7. 15. 17. 9. //

Scilab code Exa 3.3 2D Linear Convolution 1 2 3 4 5 6 7 8 9 10 11 12 13 14

// C a p t i o n : 2−D L i n e a r C o n v o l u t i o n // Example3 . 3 & Example3 . 6 & Example3 . 1 0 // p a g e 100 & p a g e 109 & p a g e 119 clc ; x =[1 ,2 ,3;4 ,5 ,6;7 ,8 ,9]; h = [3 ,4 ,5]; y = conv2d2 (x , h ) ; disp (y , ’ L i n e a r 2D c o n v o l u t i o n r e s u l t y = ’ ) // R e s u l t // L i n e a r 2D c o n v o l u t i o n r e s u l t y = // // 3. 10. 22. 22. 15. // 12. 31. 58. 49. 30. // 21. 52. 94. 76. 45.

16

Scilab code Exa 3.7 2D Linear Convolution 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

// C a p t i o n : 2−D L i n e a r C o n v o l u t i o n // Example3 . 7 // p a g e 111 clc ; x =[1 ,2;3 ,4]; h = [5 ,6;7 ,8]; y = conv2d2 (x , h ) ; disp (y , ’ L i n e a r 2D c o n v o l u t i o n r e s u l t y = ’ ) // R e s u l t // L i n e a r 2D c o n v o l u t i o n r e s u l t y = // L i n e a r 2D c o n v o l u t i o n r e s u l t y = // // 5. 16. 12. // 22. 60. 40. // 21. 52. 32

Scilab code Exa 3.8 2D Linear Convolution 1 2 3 4 5 6 7 8 9 10 11 12

// C a p t i o n : 2−D L i n e a r C o n v o l u t i o n // Example3 . 8 // p a g e 113 clc ; x =[1 ,2 ,3;4 ,5 ,6;7 ,8 ,9]; h = [1;1;1]; y = conv2d2 (x , h ) ; disp (y , ’ L i n e a r 2D c o n v o l u t i o n r e s u l t y = ’ ) // R e s u l t // L i n e a r 2D c o n v o l u t i o n r e s u l t y = // // 1 . 2. 3. // 5. 7. 9. 17

13 14 15

// // //

12. 11. 7.

15. 13. 8.

18. 15. 9.

Scilab code Exa 3.11 Linear Convolution of any signal with an impule signal given rise to the same signal 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

// C a p t i o n : L i n e a r C O n v o l u t i o n o f any s i g n a l w i t h an impulse s i g n a l gives // r i s e t o t h e same s i g n a l // Example3 . 1 1 // p a g e 121 clc ; x =[1 ,2;3 ,4]; h = 1; y = conv2d2 (x , h ) ; disp (y , ’ L i n e a r 2D c o n v o l u t i o n r e s u l t y = ’ ) // R e s u l t // L i n e a r 2D c o n v o l u t i o n r e s u l t y = // // L i n e a r 2D c o n v o l u t i o n r e s u l t y = // // 1. 2. // 3. 4.

Scilab code Exa 3.12 Circular Convolution between two 2D matrices 1 2 3 4 5 6

// C a p t i o n : C i r c u l a r C o n v o l u t i o n b e t w e e n two 2D matrices // Example3 . 1 2 // p a g e 122 clc ; x = [1 ,2;3 ,4]; h = [5 ,6;7 ,8]; 18

7 8 9 10 11 12 13 14 15 16

X = fft2d ( x ) ; // 2D FFT o f x m a t r i x H = fft2d ( h ) ; // 2D FFT o f h m a t r i x Y = X .* H ; // Element by Element m u l t i p l i c a t i o n y = ifft2d ( Y ) ; disp (y , ’ C i r c u l a r C o n v o l u t i o n R e s u l t y = ’ ) // R e s u l t // C i r c u l a r C o n v o l u t i o n R e s u l t y = // // 70. 68. // 62. 60.

Scilab code Exa 3.13 Circular Convolution exspressed as linear convolution plus alias 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

// C a p t i o n : C i r c u l a r C o n v o l u t i o n e x p r e s s e d a s l i n e a r convolution plus a l i a s // Example3 . 1 3 // p a g e 123 clc ; x = [1 ,2;3 ,4]; h = [5 ,6;7 ,8]; y = conv2d (x , h ) ; y1 = [ y (: ,1) + y (: , $ ) ,y (: ,2) ]; y2 = [ y1 (1 ,:) + y1 ($ ,:) ; y1 (2 ,:) ] disp (y , ’ L i n e a r C o n v o l u t i o n r e s u l t y= ’ ) disp ( y2 , ’ c i r c u l a r c o n v o l u t i o n e x p e s s e d a s l i n e a r convolution plus a l i a s =’) // R e s u l t // L i n e a r C o n v o l u t i o n r e s u l t y= // // 5. 16. 12. // 22. 60. 40. // 21. 52. 32. // // c i r c u l a r c o n v o l u t i o n e x p e s s e d a s l i n e a r 19

convolution plus a l i a s = 20 21 22 23

// // // //

70. 62.

68. 60.

Scilab code Exa 3.14 Linear Cross correlation of a 2D matrix 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

// C a p t i o n : l i n e a r c r o s s c o r r e l a t i o n o f a 2D m a t r i x // Example3 . 1 4 // p a g e 129 clc ; x = [3 ,1;2 ,4]; h1 = [1 ,5;2 ,3]; h2 = h1 (: , $ : -1:1) ; h = h2 ( $ : -1:1 ,:) ; y = conv2d (x , h ) disp (y , ’ L i n e a r c r o s s C o r r e l a t i o n r e s u l t y= ’ ) // R e s u l t // L i n e a r c r o s s C o r r e l a t i o n r e s u l t y= // // 9. 9. 2. // 21. 24. 9. // 10. 22. 4.

Scilab code Exa 3.15 Circular correlation between two signals 1 // C a p t i o n : C i r c u l a r 2 // Example3 . 1 5 3 // p a g e 131 4 clc ; 5 x = [1 ,5;2 ,4]; 6 h = [3 ,2;4 ,1];

c o r r e l a t i o n b e t w e e n two s i g n a l s

20

7 8 9 10 11 12 13 14 15 16 17 18

h = h (: , $ : -1:1) ; h = h ( $ : -1:1 ,:) ; X = fft2d ( x ) ; H = fft2d ( h ) ; Y = X .* H ; y = ifft2d ( Y ) ; disp (y , ’ C i r c u l a r C o r r e l a t i o n r e s u l t y= ’ ) // R e s u l t // C i r c u l a r C o r r e l a t i o n r e s u l t y= // // 37. 23. // 35. 25.

Scilab code Exa 3.16 Circular correlation between two signals 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

// C a p t i o n : C i r c u l a r c o r r e l a t i o n b e t w e e n two s i g n a l s // Example3 . 1 6 // p a g e 134 clc ; x = [5 ,10;15 ,20]; h = [3 ,6;9 ,12]; h = h (: , $ : -1:1) ; h = h ( $ : -1:1 ,:) ; X = fft2d ( x ) ; H = fft2d ( h ) ; Y = X .* H ; y = ifft2d ( Y ) ; disp (y , ’ C i r c u l a r C o r r e l a t i o n r e s u l t y= ’ ) // R e s u l t // C i r c u l a r C o r r e l a t i o n r e s u l t y= // // 300. 330. // 420. 450.

21

Scilab code Exa 3.17 Linear auto correlation of a 2D matrix 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

// C a p t i o n : l i n e a r a u t o c o r r e l a t i o n o f a 2D m a t r i x // Example3 . 1 7 // p a g e 136 clc ; x1 = [1 ,1;1 ,1]; x2 = x1 (: , $ : -1:1) ; x2 = x2 ( $ : -1:1 ,:) ; x = conv2d ( x1 , x2 ) disp (x , ’ L i n e a r a u t o C o r r e l a t i o n r e s u l t x= ’ ) // R e s u l t // L i n e a r a u t o C o r r e l a t i o n r e s u l t x= // // 1. 2. 1. // 2. 4. 2. // 1. 2. 1.

Scilab code Exa 3.18 Linear Cross correlation of a 2D matrix 1 2 3 4 5 6 7 8 9 10 11 12

// C a p t i o n : l i n e a r c r o s s c o r r e l a t i o n o f a 2D m a t r i x // Example3 . 1 8 // p a g e 141 clc ; x = [1 ,1;1 ,1]; h1 = [1 ,2;3 ,4]; h2 = h1 (: , $ : -1:1) ; h = h2 ( $ : -1:1 ,:) ; y = conv2d (x , h ) disp (y , ’ L i n e a r c r o s s C o r r e l a t i o n r e s u l t y= ’ ) // R e s u l t // L i n e a r c r o s s C o r r e l a t i o n r e s u l t y= 22

13 14 15 16

// // // //

4. 6. 2.

7. 10. 3.

3. 4. 1.

23

Chapter 4 Image Transforms

Scilab code Exa 4.4 DFT of 4x4 grayscale image 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

// C a p t i o n : 2D DFT o f 4 x4 g r a y s c a l e image // Example4 . 4 // p a g e 170 clc ; f = [1 ,1 ,1 ,1;1 ,1 ,1 ,1;1 ,1 ,1 ,1;1 ,1 ,1 ,1]; N =4; //4− p o i n t DFT kernel = dft_mtx ( N ) ; F = kernel *( f * kernel ’) ; disp (F , ’ 2D DFT o f g i v e n 2D image = ’ ) // R e s u l t // 2D DFT o f g i v e n 2D image = // // 16. 0 0 0 // 0 0 0 0 // 0 0 0 0 // 0 0 0 0

Scilab code Exa 4.5 2D DFT of 4X4 grayscale image 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

// C a p t i o n : 2D DFT o f 4 x4 g r a y s c a l e image // Example4 . 5 // p a g e 171 clc ; F = [16 ,0 ,0 ,0;0 ,0 ,0 ,0;0 ,0 ,0 ,0;0 ,0 ,0 ,0]; N =4; //4− p o i n t DFT kernel = dft_mtx ( N ) ; f = ( kernel *( F * kernel ’) ) /( N ^2) ; f = real ( f ) ; disp (f , ’ I n v e r s e 2D DFT o f t h e t r a n s f o r m e d image f = ’ ) // R e s u l t // I n v e r s e 2D DFT o f t h e t r a n s f o r m e d image f = // // 1. 1. 1. 1. // 1. 1. 1. 1. // 1. 1. 1. 1. // 1. 1. 1. 1. check Appendix AP 1 for dependency: fft2d.sce check Appendix AP 2 for dependency: ifft2d.sce

Scilab code Exa 4.6 Scilab code to intergchange phase information between two images 1 2 3 4 5

// C a p t i o n : S c i l a b c o d e t o i n t e r g c h a n g e p h a s e i n f o r m a t i o n b e t w e e n two i m a g e s // Example4 . 6 // p a g e 174 −175 clc ; close ; 25

6 a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 4 \ l e n a . png ’ ) ; 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

// SIVP t o o l b o x b = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 4 \ baboon . png ’ ) ; a = rgb2gray ( a ) ; b = rgb2gray ( b ) ; a = imresize (a ,0.5) ; b = imresize (b ,0.5) ; figure (1) ShowImage (a , ’ O r i g i n a l l e n a Image ’ ) ; // IPD t o o l b o x title ( ’ O r i g i n a l l e n a Image ’ ) ; figure (2) ShowImage (b , ’ O r i g i n a l baboon Image ’ ) ; title ( ’ O r i g i n a l baboon Image ’ ) ffta = fft2d ( double ( a ) ) ; fftb = fft2d ( double ( b ) ) ; mag_a = abs ( ffta ) ; mag_b = abs ( fftb ) ; ph_a = atan ( imag ( ffta ) , real ( ffta ) ) ; ph_b = atan ( imag ( fftb ) , real ( fftb ) ) ; newfft_a = mag_a .*( exp ( %i * ph_b ) ) ; newfft_b = mag_b .*( exp ( %i * ph_a ) ) ; rec_a = ifft2d ( newfft_a ) ; rec_b = ifft2d ( newfft_b ) ; figure (3) ShowImage ( uint8 ( rec_a ) , ’ l e n a Image a f t e r p h a s e r e v e r s a l ’ ); title ( ’ l e n a Image a f t e r p h a s e r e v e r s a l ’ ) figure (4) ShowImage ( uint8 ( rec_b ) , ’ baboon Image a f t e r p h a s e r e v e r s a l ’ ); title ( ’ baboon Image a f t e r p h a s e r e v e r s a l ’ )

Scilab code Exa 4.10 Program to compute discrete cosine transform 26

Figure 4.1: Scilab code to intergchange phase information between two images

27

1 2 3 4 5 6 7 8 9 10 11 12 13 14

// C a p t i o n : Program t o compute d i s c r e t e c o s i n e tranform // Example4 . 1 0 // p a g e 198 clc ; N =4; //DCT m a t r i x o f o r d e r f o u r X = dct_mtx ( N ) ; disp (X , ’DCT m a t r i x o f o r d e r f o u r ’ ) // R e s u l t //DCT m a t r i x o f o r d e r f o u r // // 0.5 0.5 0.5 0.5 // 0.6532815 0.2705981 − 0.2705981 − 0.6532815 // 0.5 − 0.5 − 0.5 0.5 // 0.2705981 − 0.6532815 0.6532815 − 0.2705981

Scilab code Exa 4.12 Program to perform KL tranform for the given 2D matrix 1 2 3 4 5 6 7 8 9

// C a p t i o n : Program t o p e r f o r m KL t r a n s f o r m f o r t h e g i v e n 2D m a t r i x // Example4 . 1 2 // p a g e 208 clear ; clc ; X = [4 ,3 ,5 ,6;4 ,2 ,7 ,7;5 ,5 ,6 ,7]; [m , n ]= size ( X ) ; A = []; E = []; 28

Figure 4.2: Program to compute discrete cosine transform

29

Figure 4.3: Program to compute discrete cosine transform

30

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

for i =1: n A = A + X (: , i ) ; E = E + X (: , i ) * X (: , i ) ’; end mx = A / n ; // mean m a t r i x E = E/n; C = E - mx * mx ’; // c o v a r i a n c e m a t r i x C = E [ xx ’ ] −mx∗mx ’ [V , D ] = spec ( C ) ; // e i g e n v a l u e s and e i g e n v e c t o r s d = diag ( D ) ; // d i a g o n a l e l e m e n t s od e i g e n v a l u e s [d , i ] = gsort ( d ) ; // s o r t i n g t h e e l e m e n t s o f D i n descending order for j = 1: length ( d ) T (: , j ) = V (: , i ( j ) ) ; end T =T ’ disp (d , ’ E i g e n V a l u e s a r e U = ’ ) disp (T , ’ The e i g e n v e c t o r m a t r i x T = ’ ) disp (T , ’ The KL t r a n f o r m b a s i s i s = ’ ) //KL t r a n s f o r m for i = 1: n Y (: , i ) = T * X (: , i ) ; end disp (Y , ’KL t r a n s f o r m a t i o n o f t h e i n p u t m a t r i x Y = ’ ) // R e c o n s t r u c t i o n for i = 1: n x (: , i ) = T ’* Y (: , i ) ; end disp (x , ’ R e c o n s t r u c t m a t r i x o f t h e g i v e n s a m p l e matrix X = ’ ) // R e s u l t // E i g e n V a l u e s a r e U = // 6.1963372 // 0.2147417 // 0.0264211 // The e i g e n v e c t o r m a t r i x T = // 0.4384533 0.8471005 0.3002988 // 0.4460381 − 0.4951684 0.7455591 // − 0 . 7 8 0 2 6 2 0 0.1929481 0.5949473 31

46 47 48 49 50 51 52 53 54 55 56 57

// The KL t r a n f o r m b a s i s i s = // 0.4384533 0.8471005 0.3002988 // 0.4460381 − 0.4951684 0.7455591 // − 0 . 7 8 0 2 6 2 0 0.1929481 0.5949473 // KL t r a n s f o r m a t i o n o f t h e i n p u t m a t r i x Y = // 6.6437095 4.5110551 9.9237632 10.662515 // 3.5312743 4.0755729 3.2373664 4.4289635 // 0.6254808 1.0198466 1.0190104 0.8336957 // R e c o n s t r u c t m a t r i x o f t h e g i v e n s a m p l e m a t r i x x = // 4. 3. 5. 6. // 4. 2. 7. 7. // 5. 5. 6. 7.

Scilab code Exa 4.13 Program to find the singular value decomposition of given matrix 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

// C a p t i o n : Program t o f i n d t h e s i n g u l a r v a l u e decomposition of given matrix // Example4 . 1 3 // p a g e 210 clear ; clc ; A = [1 , -2 ,3;3 ,2 , -1]; [U ,S , V ]= svd ( A ) ; A_recon = U * S *V ’; disp (U , ’U = ’ ) disp (S , ’ S = ’ ) disp (V , ’V = ’ ) disp ( A_recon , ’A m a t r i x from s v d = ’ ) // R e s u l t // U = // 32

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

// − 0 . 7 0 7 1 0 6 8 0.7071068 // 0.7071068 0.7071068 // // S = // // 4.2426407 0. // 0. 3.1622777 // // V = // // 0.3333333 0.8944272 // 0 . 6 6 6 6 6 6 7 − 2 . 7 7 6D−16 // − 0 . 6 6 6 6 6 6 7 0.4472136 // // A m a t r i x from s v d = // // 1. − 2. 3. // 3. 2. − 1.

33

0. 0.

− 0.2981424 0.7453560 0.5962848

Chapter 5 Image Enhancement

Scilab code Exa 5.5 Scilab code for brightness enhancement 1 // C a p t i o n : S c i l a b c o d e f o r b r i g h t n e s s enhancement 2 // F i g 5 . 5 3 // p a g e 246 4 clc ; 5 close ; 6 // a = i m r e a d ( ’ E : \ DIP JAYARAMAN\ C h a p t e r 5 \ p l a t e . GIF ’ ) ; 7 8 9 10 11 12 13 14 15 16

// SIVP t o o l b o x a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 4 \ baboon . png ’ ) ; a = rgb2gray ( a ) ; b = double ( a ) +50; b = uint8 ( b ) ; figure (1) ShowImage (a , ’ O r i g i n a l Image ’ ) ; title ( ’ O r i g i n a l Image ’ ) figure (2) ShowImage (b , ’ Enhanced Image ’ ) ; title ( ’ Enhanced Image ’ )

34

35

Scilab code Exa 5.7 Scilab code for brightness suppression 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

// C a p t i o n : S c i l a b c o d e f o r brightness suppression // F i g 5 . 7 // p a g e 247 clc ; close ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 4 \ baboon . png ’ ) ; a = rgb2gray ( a ) ; b = double ( a ) -50; b = uint8 ( b ) ; figure (1) ShowImage (a , ’ O r i g i n a l Image ’ ) ; title ( ’ O r i g i n a l Image ’ ) figure (2) ShowImage (b , ’ B r i g h t n e s s S u p r e s s e d Image ’ ) ; title ( ’ B r i g h t n e s s S u p r e s s e d Image ’ )

Scilab code Exa 5.9 Scilab code for Contrast Manipulation 1 2 3 4 5 6 7 8 9 10

// C a p t i o n : S c i l a b c o d e f o r Contrast Manipulation // F i g 5 . 9 // p a g e 248 clc ; close ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 4 \ l e n a . png ’ ) ; a = rgb2gray ( a ) ; b = double ( a ) *0.5; b = uint8 ( b ) c = double ( b ) *2; 36

37

11 c = uint8 ( c ) 12 figure (1) 13 ShowImage (a , ’ O r i g i n a l Image ’ ) ; 14 title ( ’ O r i g i n a l Image ’ ) 15 figure (2) 16 ShowImage (b , ’ D e c r e a s e i n C o n t r a s t ’ ) ; 17 title ( ’ D e c r e a s e i n C o n t r a s t ’ ) 18 figure (3) 19 ShowImage (c , ’ I n c r e a s e i n C o n t r a s t ’ ) ; 20 title ( ’ I n c r e a s e i n C o n t r a s t ’ )

Scilab code Exa 5.13 Scilab code to determine image negative 1 2 3 4 5 6 7 8 9 10 11 12

// C a p t i o n : S c i l a b c o d e t o d e t e r m i n e image n e g a t i v e // F i g . 5 . 1 3 // p a g e 252 clc ; close ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 5 \ l a b e l . j p g ’ ) ; k = 255 - double ( a ) ; k = uint8 ( k ) ; imshow ( a ) ; title ( ’ O r i g i n a l o n c a Image ’ ) imshow ( k ) ; title ( ’ N e g a t i v e o f O r i g i n a l Image ’ )

Scilab code Exa 5.16 Scilab code that performs threshold operation

38

Figure 5.3: Scilab code for Contrast Manipulation

Figure 5.4: Scilab code to determine image negative 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

// C a p t i o n : S c i l a b c o d e t h a t p e r f o r m s t h r e s h o l d operation // F i g 5 . 1 6 // p a g e 254 clc ; close ; a = imread ( ’E : \ D i g i t a l I m a g e P r o c e s s i n g J a y a r a m a n \ C h a p t e r 5 \ l e n a . png ’ ) ; a = rgb2gray ( a ) ; [ m n ] = size ( a ) ; t = input ( ’ E n t e r t h e t h r e s h o l d p a r a m e t e r ’ ) ; for i = 1: m for j = 1: n if ( a (i , j )
Scilab code Exa 5.20 Program performs gray level slicing without background

40

41

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

// C a p t i o n : Program p e r f o r m s g r a y l e v e l s l i c i n g without background // F i g . 5 . 2 0 // p a g e 2 5 6 clc ; x = imread ( ’E : \ D i g i t a l I m a g e P r o c e s s i n g J a y a r a m a n \ C h a p t e r 5 \ l e n a . png ’ ) ; x = rgb2gray ( x ) ; y = double ( x ) ; [m , n ]= size ( y ) ; L = max ( max ( x ) ) ; a = round ( L /2) ; b = L; for i =1: m for j =1: n if ( y (i , j ) >= a & y (i , j ) <= b ) z (i , j ) = L ; else z (i , j ) =0; end end end z = uint8 ( z ) ; figure (1) ShowImage (x , ’ O r i g i n a l Image ’ ) ; title ( ’ O r g i n a l Image ’ ) figure (2) ShowImage (z , ’ Gray L e v e l S l i c i n g ’ ) ; title ( ’ Gray L e v e l S l i c i n g w i t h o u t p r e s e r v i n g background ’ )

42

Figure 5.6: Program performs gray level slicing without background

43

Chapter 6 Image Restoration and Denoising

Scilab code Exa 6.1 Scilab code to create motion blur 1 // C a p t i o n : S c i l a b c o d e t o c r e a t e m o t i o n b l u r 2 // F i g 6 . 1 3 // p a g e 326 4 clc ; 5 close ; 6 a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 6 \humm . j p g ’ ) ; // 7 8 9 10 11 12 13 14 15 16 17 18

SIVP t o o l b o x // f i l t e r c o e f f i c i e n t s o f f s p e c i a l ( ’ motion ’ , 1 0 , 2 5 ) H =[0 ,0 ,0 ,0 ,0 ,0 ,0 ,0.0032 ,0.0449 ,0.0865 ,0.0072;... 0 ,0 ,0 ,0 ,0 ,0.0092 ,0.0509 ,0.0925 ,0.0629 ,0.0213 ,0;... 0 ,0 ,0 ,0.0152 ,0.0569 ,0.0985 ,0.0569 ,0.0152 ,0 ,0 ,0;... 0 ,0.0213 ,0.0629 ,0.0925 ,0.0509 ,0.0092 ,0 ,0 ,0 ,0 ,0;... 0.0072 ,0.0865 ,0.0449 ,0.0032 ,0 ,0 ,0 ,0 ,0 ,0 ,0]; Motion_Blur = imfilter (a , H ) ; Motion_Blur = uint8 ( Motion_Blur ) ; ShowImage (a , ’ o r i g i n a l Image ’ ) title ( ’ o r i g i n a l Image ’ ) figure ShowImage ( Motion_Blur , ’ Motion B l u r r e d Image ’ ) 44

Figure 6.1: Scilab code to create motion blur 19

title ( ’ 10 x25 Motion B l u r r e d Image ’ )

check Appendix AP 1 for dependency: fft2d.sce check Appendix AP 2 for dependency: ifft2d.sce

Scilab code Exa 6.5 Scilab code performs inverse filtering 1 2 3 4 5 6 7

// C a p t i o n : S c i l a b c o d e p e r f o r m s i n v e r s e f i l t e r i n g // Degrade t h e image by means o f a known b l u r // Apply i n v e r s e f i l t e r t o t h e b l u r r e d image and s e e t h e r e s t o r e d image // F i g 6 . 5 // p a g e 330 clc ; close ; 45

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37

x = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 6 \ f l o w e r 2 . j p g ’ ) ; x = double ( rgb2gray ( x ) ) ; [ M N ]= size ( x ) ; h = zeros (M , N ) ; for i = 1:11 for j = 1:11 h (i , j ) = 1/121; end end sigma = sqrt (4*10^( -7) ) ; freqx = fft2d ( x ) ; // F o u r i e r t r a n s f o r m o f i n p u t image freqh = fft2d ( h ) ; // F o u r i e r t r a n s f o r m o f d e g r a d a t i o n y = real ( ifft2d ( freqh .* freqx ) ) ; freqy = fft2d ( y ) ; powfreqx = freqx .^2/( M * N ) ; alpha = 0.5; // I n d i c a t e s i n v e r s e f i l t e r freqg = (( freqh . ’) ’) .* abs ( powfreqx ) ./( abs ( freqh .^2) .* abs ( powfreqx ) + alpha * sigma ^2) ; Resfreqx = freqg .* freqy ; Resa = real ( ifft2d ( Resfreqx ) ) ; x = uint8 ( x ) ; y = uint8 ( y ) ; Resa = uint8 ( Resa ) ShowImage (x , ’ O r i g i n a l Image ’ ) title ( ’ O r i g i n a l Image ’ ) figure ShowImage (y , ’ Degraded Image ’ ) title ( ’ Degraded Image ’ ) figure ShowImage ( Resa , ’ R e s t o r e d Image ’ ) title ( ’ R e s t o r e d Image ’ )

check Appendix AP 1 for dependency: fft2d.sce check Appendix AP 2 for dependency: 46

47

ifft2d.sce

Scilab code Exa 6.7 Scilab code performs inverse filtering 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

// C a p t i o n : S c i l a b c o d e p e r f o r m s i n v e r s e f i l t e r i n g // Degrade t h e image by means o f a known b l u r and white noise // The image i s d e g r a d e d a s w e l l a s c o r r u p t e d by noise // Apply i n v e r s e f i l t e r t o r e s t o r e t h e image // F i g 6 . 7 // p a g e 332 clc ; close ; x = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 6 \ f l o w e r 2 . j p g ’ ) ; x = double ( rgb2gray ( x ) ) ; [ M N ]= size ( x ) ; h = zeros (M , N ) ; for i = 1:11 for j = 1:11 h (i , j ) = 1/121; end end sigma = sqrt (4*10^( -7) ) ; freqx = fft2d ( x ) ; // F o u r i e r t r a n s f o r m o f i n p u t image freqh = fft2d ( h ) ; // F o u r i e r t r a n s f o r m o f d e g r a d a t i o n y = real ( ifft2d ( freqh .* freqx ) ) +10* rand (M ,N , ’ n o r m a l ’ ) ; freqy = fft2d ( y ) ; powfreqx = freqx .^2/( M * N ) ; alpha = 0.5; // I n d i c a t e s i n v e r s e f i l t e r freqg = (( freqh . ’) ’) .* abs ( powfreqx ) ./( abs ( freqh .^2) .* abs ( powfreqx ) + alpha * sigma ^2) ; Resfreqx = freqg .* freqy ; 48

49

27 Resa = real ( ifft2d ( Resfreqx ) ) ; 28 x = uint8 ( x ) ; 29 y = uint8 ( y ) ; 30 Resa = uint8 ( Resa ) 31 ShowImage (x , ’ O r i g i n a l Image ’ ) 32 title ( ’ O r i g i n a l Image ’ ) 33 figure 34 ShowImage (y , ’ Degraded+n o i s e Image ’ ) 35 title ( ’ Degraded+n o i s e Image ’ ) 36 figure 37 ShowImage ( Resa , ’ R e s t o r e d Image ’ ) 38 title ( ’ R e s t o r e d Image ’ )

check Appendix AP 1 for dependency: fft2d.sce check Appendix AP 2 for dependency: ifft2d.sce

Scilab code Exa 6.9 Scilab code performs Pseudo inverse filtering 1 2 3 4 5 6 7 8 9

// C a p t i o n : S c i l a b c o d e p e r f o r m s Pseudo i n v e r s e filtering // Degrade t h e image by means o f a known b l u r and white noise // The image i s d e g r a d e d a s w e l l a s c o r r u p t e d by noise // Apply Pseudo i n v e r s e f i l t e r t o r e s t o r e t h e image // F i g 6 . 9 // p a g e 333 clc ; close ; x = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 6 \ f l o w e r 2 . j p g ’ ) ; 50

51

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

x = double ( rgb2gray ( x ) ) ; [ M N ]= size ( x ) ; h = zeros (M , N ) ; for i = 1:11 for j = 1:11 h (i , j ) = 1/121; end end mask_b = ones (11 ,11) /121; [ m1 , n1 ] = size ( mask_b ) ; Thr_Freq = 0.2; freqx = fft2d ( x ) ; // F o u r i e r t r a n s f o r m o f i n p u t image freqh = fft2d ( h ) ; // F o u r i e r t r a n s f o r m o f d e g r a d a t i o n y = real ( ifft2d ( freqh .* freqx ) ) +25* rand (M ,N , ’ n o r m a l ’ ) ; freqy = fft2d ( y ) ; psf = zeros (M , N ) ; psf ( M /2+1 -( m1 -1) /2: M /2+1+( m1 -1) /2 , N /2+1 -( n1 -1) /2: N /2+1+( n1 -1) /2) = mask_b ; psf = fftshift ( psf ) ; freq_res = fft2d ( psf ) ; Inv_filt = freq_res ./(( abs ( freq_res ) ) .^2+ Thr_Freq ) ; z = real ( ifft2d ( freqy .* Inv_filt ) ) ; x = uint8 ( x ) ; y = uint8 ( y ) ; z = uint8 ( z ) ShowImage (x , ’ O r i g i n a l Image ’ ) title ( ’ O r i g i n a l Image ’ ) figure ShowImage (y , ’ Degraded+n o i s e Image ’ ) title ( ’ Degraded+n o i s e Image ’ ) figure ShowImage (z , ’ R e s t o r e d Image ’ ) title ( ’ R e s t o r e d Image ’ )

check Appendix AP 1 for dependency: 52

53

fft2d.sce check Appendix AP 2 for dependency: ifft2d.sce

Scilab code Exa 6.13 Scilab code to perform wiener filtering of the corrupted image 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

// C a p t i o n : S c i l a b c o d e t o p e r f o r m w i e n e r f i l t e r i n g o f t h e c o r r u p t e d image // F i g 6 . 1 3 // Page 339 close ; clc ; x = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 6 \ f l o w e r 2 . j p g ’ ) ; // SIVP t o o l b o x x = double ( rgb2gray ( x ) ) ; sigma = 50; Gamma = 1; alpha = 1; // I t i n d i c a t e s Wiener f i l t e r [ M N ]= size ( x ) ; h = zeros (M , N ) ; for i = 1:5 for j = 1:5 h (i , j ) = 1/25; end end Freqa = fft2d ( x ) ; Freqh = fft2d ( h ) ; y = real ( ifft2d ( Freqh .* Freqa ) ) // image d e g r a d a t i o n y = y +25* rand (M ,N , ” n o r m a l ” ) ; // Adding random n o i s e with normal d i s t r i b u t i o n Freqy = fft2d ( y ) ; Powy = abs ( Freqy ) .^2/( M * N ) ; sFreqh = Freqh .*( abs ( Freqh ) >0) +1/ Gamma *( abs ( Freqh ) ==0) ; 54

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

iFreqh = 1/ sFreqh ; iFreqh = iFreqh ’.*( abs ( Freqh ) * Gamma >1) + Gamma * abs ( sFreqh ) * iFreqh *( abs ( sFreqh ) * Gamma <=1) ; iFreqh = iFreqh /( max ( max ( abs ( iFreqh ) ) ) ) ; Powy = Powy .*( Powy > sigma ^2) + sigma ^2*( Powy <= sigma ^2) ; Freqg = iFreqh .*( Powy - sigma ^2) ./( Powy -(1 - alpha ) * sigma ^2) ; ResFreqa = Freqg .* Freqy ; Resa = real ( ifft2d ( ResFreqa ) ) ; x = uint8 ( x ) ; y = uint8 ( y ) ; Resa = uint8 ( Resa ) ; ShowImage (x , ’ O r i g i n a l Image ’ ) title ( ’ O r i g i n a l Image ’ ) figure ShowImage (y , ’ Degraded Image ’ ) title ( ’ Degraded Image ’ ) figure ShowImage ( Resa , ’ R e s t o r e d Image ’ ) title ( ’ R e s t o r e d Image ’ )

Scilab code Exa 6.18 Scilab code to Perform Average Filtering operation 1 2 3 4 5 6 7

// C a p t i o n : S c i l a b c o d e t o P e r f o r m A v e r a g e F i l t e r i n g operation // F i g 6 . 1 8 // p a g e 349 clc ; close ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 6 \ l e n n a . j p g ’ ) ; // SIVP t o o l b o x a = imnoise (a , ’ s a l t & p e p p e r ’ , 0.2) ; //Add s a l t &p e p p e r n o i s e t o t h e image 55

56

8 a = double ( a ) ; 9 [ m n ]= size ( a ) ; 10 N = input ( ’ e n t e r t h e window s i z e = ’ ) ; // The window s i z e 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

can be 3 x3 , 5 x 5 e t c Start =( N +1) /2; Out_Imag = a ; for i = Start :( m - Start +1) for j = Start :( n - Start +1) limit =( N -1) /2; Sum =0; for k = - limit : limit , for l = - limit : limit , Sum = Sum + a ( i +k , j + l ) ; end end Out_Imag (i , j ) = Sum /( N * N ) ; end end a = uint8 ( a ) ; Out_Imag = uint8 ( Out_Imag ) ; ShowImage (a , ’ o r i g i n a l Image ’ ) title ( ’ N o i s y Image ’ ) figure ShowImage ( Out_Imag , ’ a v e r a g e f i l t e r e d Image ’ ) title ( ’ 5 x5 a v e r a g e f i l t e r e d Image ’ ) ;

Scilab code Exa 6.21 Scilab code to Perform median filtering 1 // C a p t i o n : S c i l a b c o d e t o P e r f o r m median 2 // F i g 6 . 2 1 3 // p a g e 352 4 clc ; 5 close ;

57

filtering

Figure 6.7: Scilab code to Perform Average Filtering operation 6 c = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 6 \ cameraman . j p g ’ 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

) ; // SIVP t o o l b o x N = input ( ’ E n t e r t h e window s i z e ’ ) ; a = double ( imnoise (c , ’ s a l t & p e p p e r ’ ,0.2) ) ; [m , n ] = size ( a ) ; b = a; if ( modulo (N ,2) ==1) Start = ( N +1) /2; End = Start ; limit1 = (N -1) /2; limit2 = limit1 ; else Start = N /2; End = Start +1; limit1 = ( N /2) -1; limit2 = limit1 +1; end for i = Start :( m - End +1) for j = Start :( n - End +1) I =1; for k = - limit1 : limit2 for l = - limit1 : limit2 mat ( I ) = a ( i +k , j +1) I = I +1; end 58

30 end 31 mat = gsort ( mat ) ; 32 if ( modulo (N ,2) ==1) 33 b (i , j ) = ( mat ((( N ^2) +1) /2) ) ; 34 else 35 b (i , j ) = ( mat (( N ^2) /2) + mat ((( N ^2) /2) +1) ) /2; 36 end 37 end 38 end 39 a = uint8 ( a ) ; 40 b = uint8 ( b ) ; 41 figure 42 ShowImage (c , ’ O r i g i n a l Image ’ ) 43 title ( ’ O r i g i n a l Image ’ ) 44 figure 45 ShowImage (a , ’ n o i s y image ’ ) 46 title ( ’ n o i s y image ’ ) 47 figure 48 ShowImage (b , ’ Median F i l t e r e d Image ’ ) 49 title ( ’ 5 x5 Median F i l t e r e d Image ’ )

check Appendix AP 3 for dependency: Func_medianall.sci

Scilab code Exa 6.23 Scilab code to Perform median filtering of colour image 1 2 3 4 5

// C a p t i o n : S c i l a b c o d e t o P e r f o r m median f i l t e r i n g o f c o l o u r image // F i g 6 . 2 3 ( a ) // p a g e 353 clc ; close ; 59

Figure 6.8: Scilab code to Perform median filtering 6 a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 6 \ p e p p e r s . png ’ ) ; 7 8 9 10 11 12 13 14 15 16 17 18 19

// SIVP t o o l b o x N = input ( ’ e n t e r t h e window s i z e ’ ) ; b = imresize (a ,[256 ,256]) ; b = imnoise (b , ’ s a l t & p e p p e r ’ ,.1) ; [ m n ]= size ( b ) ; R = b (: ,: ,1) ; G = b (: ,: ,2) ; B = b (: ,: ,3) ; Out_R = Func_medianall (R , N ) ; // A p p l y i n g Median f i l t e r to R plane Out_G = Func_medianall (G , N ) ; // A p p l y i n g Median f i l t e r to G plane Out_B = Func_medianall (B , N ) ; // A p p l y i n g Median f i l t e r to B plane Out_Image (: ,: ,1) = Out_R ; Out_Image (: ,: ,2) = Out_G ; Out_Image (: ,: ,3) = Out_B ; 60

Figure 6.9: Scilab code to Perform median filtering of colour image 20 b = uint8 ( b ) ; 21 Out_Image = uint8 ( Out_Image ) ; 22 // ShowColorImage ( b , ’ n o i s e added ’ ) 23 // t i t l e ( ’ n o i s e added ’ ) 24 figure 25 ShowColorImage ( Out_Image , ’ 3 x3 median 26 title ( ’ 3 x3 median f i l t e r e d ’ )

filtered ’)

Scilab code Exa 6.24 Scilab code to Perform Trimmed Average Filter 1 2 3 4 5 6

// C a p t i o n : S c i l a b c o d e t o P e r f o r m Trimmed A v e r a g e Filter // Alpha trimmed a v e r a g e f i l t e r // F i g 6 . 2 4 // p a g e 355 clc ; close ; 61

7 c = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 6 \ l e n n a . j p g ’ ) ; //

SIVP t o o l b o x 8 s = 1; // s d e n o t e s t h e number o f 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

v a l u e s t o be l e f t

i n t h e end r = 1; N = 9; // 3 x3 window a = double ( imnoise (c , ’ g a u s s i a n ’ ) ) ; [m , n ] = size ( a ) ; b = zeros (m , n ) ; for i = 2: m -1 for j = 2: n -1 mat = [ a (i , j ) ,a (i ,j -1) ,a (i , j +1) ,a (i -1 , j ) ,a ( i +1 , j ) ,a (i -1 ,j -1) ,... a (i -1 , j +1) ,a (i -1 , j +1) ,a ( i +1 , j +1) ]; sorted_mat = gsort ( mat ) ; Sum =0; for k = r + s :( N - s ) Sum = Sum + mat ( k ) ; end b (i , j ) = Sum /( N -r - s ) ; end end a = uint8 ( a ) ; b = uint8 ( b ) ; // f i g u r e // imshow ( c ) // t i t l e ( ’ O r i g i n a l Image ’ ) figure ShowImage (a , ’ n o i s y image ’ ) title ( ’ n o i s y image ’ ) figure ShowImage (b , ’ Trimmed A v e r a g e F i l t e r e d Image ’ ) title ( ’ Trimmed A v e r a g e F i l t e r e d Image ’ )

62

Figure 6.10: Scilab code to Perform Trimmed Average Filter

63

Chapter 7 Image Segmentation

Scilab code Exa 7.23 Scilab code for Differentiation of Gaussian function 1 2 3 4 5 6 7 8 9 10 11 12 13

// C a p t i o n : S c i l a b c o d e f o r D i f f e r e n t i a t i o n o f Gaussian f u n c t i o n // F i g 7 . 2 3 // p a g e 3 8 8 clc ; close ; sigma = input ( ’ E n t e r t h e v a l u e o f s i g m a : ’ ) i = -10:.1:10; j = -10:.1:10; r = sqrt ( i .* i + j .* j ) ; y =(1/( sigma ^2) ) *((( r .* r ) / sigma ^2) -1) .* exp ( - r .* r /2* sigma ^2) ; plot (i , y ) legend ( sprintf ( ’ The s i g m a v a l u e i s %g ’ , sigma ) ) xtitle ( ’ D i f f e r e n t i a t i o n o f G a u s s i a n f u n c t i o n ’ )

64

Figure 7.1: Scilab code for Differentiation of Gaussian function

65

Figure 7.2: Scilab code for Differentiation of Gaussian function

66

Figure 7.3: Scilab code for Differentiation of Gaussian Filter function

Scilab code Exa 7.25 Scilab code for Differentiation of Gaussian Filter function 1 2 3 4 5 6 7 8 9

// C a p t i o n : S c i l a b c o d e f o r D i f f e r e n t i a t i o n o f Gaussian F i l t e r f u n c t i o n // F i g 7 . 2 5 // p a g e 3 8 9 clc ; close ; sigma1 = input ( ’ E n t e r t h e v a l u e o f s i g m a 1 : ’ ) sigma2 = input ( ’ E n t e r t h e v a l u e o f s i g m a 2 : ’ ) i = -10:.1:10; j = -10:.1:10; 67

Figure 7.4: Scilab code for Differentiation of Gaussian Filter function 10 r = sqrt ( i .* i + j .* j ) ; 11 y1 = (1/( sigma1 ^2) ) *((( r .* r ) / sigma1 ^2) -1) .* exp ( - r .* r 12 13 14 15 16 17 18 19

/2* sigma1 ^2) ; y2 = (1/( sigma2 ^2) ) *((( r .* r ) / sigma2 ^2) -1) .* exp ( - r .* r /2* sigma2 ^2) ; y = y1 - y2 ; plot (i , y ) xtitle ( ’ Shape o f DOG F i l t e r ’ ) // R e s u l t // E n t e r t h e v a l u e o f s i g m a 1 : 4 // E n t e r t h e v a l u e o f s i g m a 2 : 1 //

68

Scilab code Exa 7.27 Scilab code for Edge Detection using Different Edge detectors 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

// C a p t i o n : S c i l a b c o d e f o r Edge D e t e c t i o n u s i n g D i f f e r e n t Edge d e t e c t o r s // [ 1 ] . S o b e l [ 2 ] . P r e w i t t [ 3 ] . Log [ 4 ] . Canny // F i g 7 . 2 7 // p a g e 3 8 9 close ; clc ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 7 \ s a i l i n g . j p g ’ ) ; a = rgb2gray ( a ) ; c = edge (a , ’ s o b e l ’ ) ; d = edge (a , ’ p r e w i t t ’ ) ; e = edge (a , ’ l o g ’ ) ; f = edge (a , ’ canny ’ ) ; ShowImage (a , ’ O r i g i n a l Image ’ ) title ( ’ O r i g i n a l Image ’ ) figure ShowImage (c , ’ S o b e l ’ ) title ( ’ S o b e l ’ ) figure ShowImage (d , ’ P r e w i t t ’ ) title ( ’ P r e w i t t ’ ) figure ShowImage (e , ’ Log ’ ) title ( ’ Log ’ ) figure ShowImage (f , ’ Canny ’ ) title ( ’ Canny ’ )

69

Figure 7.5: Scilab code for Edge Detection using Different Edge detectors

70

Scilab code Exa 7.30 Scilab code to perform watershed transform 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

// C a p t i o n : S c i l a b c o d e t o p e r f o r m w a t e r s h e d transform // F i g 7 . 3 0 // Page396 clc ; close ; b = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 7 \ t e a s e t . png ’ ) ; a = rgb2gray ( b ) ; global EDGE_SOBEL ; Gradient = EdgeFilter (a , EDGE_SOBEL ) ; Threshold1 = CalculateOtsuThreshold ( Gradient ) ; // determine a threshold EdgeImage = ~ SegmentByThreshold ( Gradient , Threshold1 ) ; DistanceImage = DistanceTransform ( EdgeImage ) ; Threshold2 = CalculateOtsuThreshold ( DistanceImage ) // d e t e r m i n e a t h r e s h o l d ThresholdImage = SegmentByThreshold ( DistanceImage , Threshold2 ) ; MarkerImage = SearchBlobs ( ThresholdImage ) ; SegmentedImage = Watershed ( Gradient , MarkerImage ) ; figure ShowColorImage (b , ’ t e a s e t ’ ) title ( ’ t e a s e t . png ’ ) figure ColorMapLength = length ( unique ( SegmentedImage ) ) ; ShowImage ( SegmentedImage , ’ R e s u l t o f Watershed T r a n s f o r m ’ , jetcolormap ( ColorMapLength ) ) ;

71

Figure 7.6: Scilab code to perform watershed transform

72

Chapter 8 Object Recognition

Scilab code Exa 8.4 To verify the given matrix is a covaraince matrix 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

// C a p t i o n : To v e r i f y t h e g i v e n m a t r i x i s a covaraince matrix // Problem 4 // p a g e 4 3 8 close ; clear all ; clc ; K = [37 , -15; -15 ,37]; evals = spec ( K ) ; evals = gsort ( evals ) ; disp ( evals , ’ E i g e n V a l u e s a r e = ’ ) if ( evals == abs ( evals ) ) then disp ( ’ Both t h e e i g e n v a l u e s a r e non−n e g a t i v e and the g i v e n matrix i s a c o v a r i a n c e matrix ’ ); else disp ( ’ non−c o v a r i a n c e m a t r i x ’ ) end

Scilab code Exa 8.5 To compute the covariance of the given 2D data 73

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

// C a p t i o n : To compute t h e c o v a r i a n c e o f t h e g i v e n 2D data // Problem 5 // p a g e 4 3 9 close ; clear all ; clc ; X1 = [2 ,1] ’; X2 = [3 ,2] ’; X3 = [2 ,3] ’; X4 = [1 ,2] ’; X = [ X1 , X2 , X3 , X4 ]; disp (X , ’X= ’ ) ; [M , N ] = size ( X ) ; //M=rows , N = c o l u m n s for i =1: N m ( i ) = mean ( X (: , i ) ) ; A (: , i ) = X (: , i ) -m ( i ) ; end m = m ’; disp (m , ’ mean = ’ ) ; K = A ’* A ; K = K /( M -1) ; disp (K , ’ The C o v a r a i n c e m a t i x i s K = ’ ) // R e s u l t //X= // 2. 3. 2. 1. // 1. 2. 3. 2. // mean = // 1.5 2.5 2.5 1.5 // // The C o v a r a i n c e m a t i x i s K = // 0.5 0.5 − 0.5 − 0.5 // 0.5 0.5 − 0.5 − 0.5 // − 0 . 5 − 0 . 5 0.5 0.5 // − 0 . 5 − 0 . 5 0.5 0.5

74

Scilab code Exa 8.9 Develop a perceptron AND function with bipolar inputs and targets 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

// C a p t i o n : D e v e l o p a p e r c e p t r o n AND f u n c t i o n w i t h b i p o l a r i n p u t s and t a r g e t s // Problem 9 // p a g e 4 4 1 close ; clear all ; clc ; X1 = [1 , -1 ,1 , -1]; //X1 and X2 a r e i n p u t v e c t o r s t o AND f u n c t i o n X2 = [1 ,1 , -1 , -1]; // b = [ 1 , 1 , 1 , 1 ] ; // B i a s i n g v e c t o r T = [1 , -1 , -1 , -1]; // T a r g e t v e c t o r f o r AND f u n c t i o n W1 = 0; // W e i g h t s a r e i n i t i a l i z e d W2 = 0; b = 0; // b i a s i n i t i a l i z e d alpha = 1; // l e a r n i n g r a t e for i = 1: length ( X1 ) Yin ( i ) = b + X1 ( i ) * W1 + X2 ( i ) * W2 ; if ( Yin ( i ) >=1) Y ( i ) =1; elseif (( Yin ( i ) <1) &( Yin ( i ) >= -1) ) Y ( i ) =0; elseif ( Yin ( i ) < -1) Y ( i ) = -1; end disp ( Yin ( i ) , ’ Yin= ’ ) disp ( Y ( i ) , ’Y= ’ ) if ( Y ( i ) ~= T ( i ) ) b = b + alpha * T ( i ) ; W1 = W1 + alpha * T ( i ) * X1 ( i ) ; W2 = W2 + alpha * T ( i ) * X2 ( i ) ; 75

30 disp (b , ’ b= ’ ) 31 disp ( W1 , ’W1= ’ ) 32 disp ( W2 , ’W2= ’ ) 33 end 34 end 35 disp ( ’ F i n a l W e i g h t s a f t e r one i t e r a t i o n 36 disp (b , ’ B i a s Weigth b= ’ ) 37 disp ( W1 , ’W1= ’ ) 38 disp ( W2 , ’W2= ’ )

76

are ’)

Chapter 9 Image Compression

Scilab code Exa 9.9 Program performs Block Truncation Coding BTC 1 2 3 4 5 6 7

8 9 10 11 12 13 14 15 16 17

// C a p t i o n : Program p e r f o r m s B l o c k T r u n c a t i o n Coding ( BTC) // Example 9 . 9 // p a g e 5 1 2 close ; clear all ; clc ; x = [65 ,75 ,80 ,70;72 ,75 ,82 ,68;84 ,72 ,62 ,65;66 ,68 ,72 ,80]; disp (x , ’ O r i g i n a l B l o c k i s x = ’ ) [ m1 n1 ]= size ( x ) ; blk = input ( ’ E n t e r t h e b l o c k s i z e : ’ ) ; for i = 1 : blk : m1 for j = 1 : blk : n1 y = x ( i : i +( blk -1) ,j : j +( blk -1) ) ; m = mean ( mean ( y ) ) ; disp (m , ’ mean v a l u e i s m = ’ ) sig = std2 ( y ) ; disp ( sig , ’ S t a n d a r d d e v i a t i o n o f t h e b l o c k i s =’) 77

18 19 20 21 22 23 24 25 26 27

b = y > m ; // t h e b i n a r y b l o c k disp (b , ’ B i n a r y a l l o c a t i o n m a t r i x i s B= ’ ) K = sum ( sum ( b ) ) ; disp (K , ’ number o f o n e s = ’ ) if ( K ~= blk ^2 ) & ( K ~= 0) ml = m - sig * sqrt ( K /(( blk ^2) -K ) ) ; disp ( ml , ’ The v a l u e o f a = ’ ) mu = m + sig * sqrt ((( blk ^2) -K ) / K ) ; disp ( mu , ’ The v a l u e o f b = ’ ) x ( i : i +( blk -1) , j : j +( blk -1) ) = b * mu +(1 - b ) * ml ; end

28 29 end 30 end 31 disp ( round ( x ) , ’ R e c o n s t r u c t e d B l o c k i s x = ’ ) 32 // R e s u l t 33 // O r i g i n a l B l o c k i s x = 34 // 35 // 65. 75. 80. 70. 36 // 72. 75. 82. 68. 37 // 84. 72. 62. 65. 38 // 66. 68. 72. 80. 39 // 40 // E n t e r t h e b l o c k s i z e : 4 41 // mean v a l u e i s m = 7 2 . 2 5 42 // S t a n d a r d d e v i a t i o n o f t h e b l o c k i s = 6 . 6 2 8 2 2 2 5 43 // B i n a r y a l l o c a t i o n m a t r i x i s B= 44 // 45 // F T T F 46 // F T T F 47 // T F F F 48 // F F F T 49 // 50 // number o f o n e s = 6 51 // The v a l u e o f a = 6 7 . 1 1 5 8 0 1 52 // The v a l u e o f b = 8 0 . 8 0 6 9 9 8 53 // R e c o n s t r u c t e d B l o c k i s x = 54 //

78

55 56 57 58

// // // //

67. 67. 81. 67.

81. 81. 67. 67.

81. 81. 67. 67.

67. 67. 67. 81.

Scilab code Exa 9.59 Program performs Block Truncation Coding 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

// C a p t i o n : Program p e r f o r m s B l o c k T r u n c a t i o n Coding ( BTC) by c h o o s i n g d i f f e r e n t // b l o c k s i z e s // F i g . 9 . 5 9 : MATLAB Example1 // p a g e 5 1 4 close ; clc ; x = imread ( ’E : \ D i g i t a l I m a g e P r o c e s s i n g J a y a r a m a n \ C h a p t e r 9 \ l e n n a . j p g ’ ) ; // SIVP t o o l b o x // x=i m r e s i z e ( x , [ 2 5 6 2 5 6 ] ) ; x1 = x ; x = double ( x ) ; [ m1 n1 ]= size ( x ) ; blk = input ( ’ E n t e r t h e b l o c k s i z e : ’ ) ; for i = 1 : blk : m1 for j = 1 : blk : n1 y = x ( i : i +( blk -1) ,j : j +( blk -1) ) ; m = mean ( mean ( y ) ) ; sig = std2 ( y ) ; b = y > m ; // t h e b i n a r y b l o c k K = sum ( sum ( b ) ) ; if ( K ~= blk ^2 ) & ( K ~= 0) ml = m - sig * sqrt ( K /(( blk ^2) -K ) ) ; mu = m + sig * sqrt ((( blk ^2) -K ) / K ) ; x ( i : i +( blk -1) , j : j +( blk -1) ) = b * mu +(1 - b ) * ml ; end end 79

Figure 9.1: Program performs Block Truncation Coding 26 end 27 // imshow ( u i n t 8 ( x ) ) 28 // t i t l e ( ’ R e c o n s t r u c t e d Image ’ ) 29 x = uint8 ( x ) ; 30 figure (1) 31 imshow ( x1 ) 32 title ( ’ O r i g i n a l Image ’ ) ; // IPD t o o l b o x 33 figure (2) 34 ShowImage (x , ’ R e c o n s t r u c t e d Image ’ ) ; // IPD t o o l b o x 35 title ( ’ B l o c k S i z e = 8 ’ )

80

Chapter 10 Binary Image Processing

Scilab code Exa 10.17 Scilab Code for dilation and erosion process 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

// C a p t i o n : S c i l a b Code f o r d i l a t i o n and e r o s i o n process // F i g . 1 0 . 1 7 // Page553 close ; clear all ; clc ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 0 \ morph1 . bmp ’ ) ; // SIVP t o o l b o x // b = [ 1 , 1 , 1 ; 1 , 1 , 1 ; 1 , 1 , 1 ] ; StructureElement = CreateStructureElement ( ’ s q u a r e ’ , 3) ; a1 = DilateImage (a , StructureElement ) ; a2 = ErodeImage (a , StructureElement ) ; // D i s p l a y i n g o r i g i n a l Image // imshow ( a ) figure (1) ShowImage (a , ’ O r i g i n a l Image ’ ) ; // D i s p l a y i n g D i l a t e d Image // imshow ( a1 ) figure (2) 81

Figure 10.1: Scilab Code for dilation and erosion process 19 20 21 22 23 24 25

ShowImage ( a1 , ’ D i l a t e d Image ’ ) ; xtitle ( ’ D i l a t e d Image ’ ) // D i s p l a y i n g Eroded Image // imshow ( a2 ) figure (3) ShowImage ( a2 , ’ Eroded Image ’ ) ; xtitle ( ’ Eroded Image ’ )

Scilab code Exa 10.19 Scilab Code to perform an opening and closing operation on the image 1 2 3 4 5 6 7

// C a p t i o n : S c i l a b Code t o p e r f o r m an o p e n i n g and c l o s i n g o p e r a t i o n on t h e image // F i g . 1 0 . 1 9 // Page555 close ; clear all ; clc ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 0 \ morph2 . bmp ’ ) ; 82

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

// SIVP t o o l b o x // b = [ 1 , 1 , 1 ; 1 , 1 , 1 ; 1 , 1 , 1 ] ; StructureElement = CreateStructureElement ( ’ s q u a r e ’ , 3) ; // Opening i s done by f i r s t a p p l y i n g e r o s i o n and t h e n d i l a t i o n o p e r a t i o n s on image b1 = ErodeImage (a , StructureElement ) ; b2 = DilateImage ( b1 , StructureElement ) ; // C l o s i n g i s done by f i r s t a p p l y i n g d i l a t i o n and t h e n e r o s i o n o p e r a t i o n on image a1 = DilateImage (a , StructureElement ) ; a2 = ErodeImage ( a1 , StructureElement ) ; // D i s p l a y i n g o r i g i n a l Image figure (1) ShowImage (a , ’ O r i g i n a l Image ’ ) ; // D i s p l a y i n g Opened Image figure (2) ShowImage ( b2 , ’ Opened Image ’ ) ; xtitle ( ’ Opened Image ’ ) // D i s p l a y i n g C l o s e d Image figure (3) ShowImage ( a2 , ’ C l o s e d Image ’ ) ; xtitle ( ’ C l o s e d Image ’ )

83

Figure 10.2: Scilab Code to perform an opening and closing operation on the image

84

Chapter 11 Colur Image Processing

Scilab code Exa 11.4 Read an RGB image and extract the three colour components red green blue 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

// C a p t i o n : Read an RGB image and e x t r a c t t h e t h r e e c o l o u r c o m p o n e n ts : red , g r e e n // and b l u e // F i g . 1 1 . 4 : MATLAB Example1 // p a g e 5 8 8 clc ; close ; RGB = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 1 \ p e p p e r s . png ’ ) ; // SIVP t o o l b o x R = RGB ; G = RGB ; B = RGB ; R (: ,: ,2) =0; R (: ,: ,3) =0; G (: ,: ,1) =0; G (: ,: ,3) =0; B (: ,: ,1) =0; B (: ,: ,2) =0; figure (1) ShowColorImage ( RGB , ’ O r i g i n a l C o l o r Image ’ ) ; // IPD 85

19 20 21 22 23 24 25

toolbox title ( ’ O r i g i n a l C o l o r Image ’ ) ; figure (2) ShowColorImage (R , ’ Red Component ’ ) ; figure (3) ShowColorImage (G , ’ Green Component ’ ) ; figure (4) ShowColorImage (B , ’ B l u e Component ’ ) ;

Scilab code Exa 11.12 Read a Colour image and separate the colour image into red green and blue planes 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

// C a p t i o n : Read a C o l o u r image and s e p a r a t e t h e c o l o u r image i n t o : red , g r e e n // and b l u e p l a n e s // F i g . 1 1 . 1 2 : MATLAB Example2 // p a g e 5 9 2 clc ; close ; RGB = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 1 \ p e p p e r s . png ’ ) ; // SIVP t o o l b o x a1 = RGB ; b1 = RGB ; c1 = RGB ; a1 (: ,: ,1) =0; b1 (: ,: ,2) =0; c1 (: ,: ,3) =0; figure (1) ShowColorImage ( RGB , ’ O r i g i n a l C o l o r Image ’ ) ; // IPD toolbox figure (2) ShowColorImage ( a1 , ’ Red M i s s i n g ’ ) ; figure (3) 86

Figure 11.1: Read an RGB image and extract the three colour components red green blue

87

19 20 21

ShowColorImage ( b1 , ’ Green M i s s i n g ’ ) ; figure (4) ShowColorImage ( c1 , ’ B l u e M i s s i n g ’ ) ;

Scilab code Exa 11.16 Compute the histogram of the colour image 1 // C a p t i o n : Compute t h e h i s t o g r a m o f t h e c o l o u r image 2 // F i g . 1 1 . 1 6 : MATLAB Example3 3 // p a g e 5 9 5 4 clc ; 5 close ; 6 I = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 1 \ l a v e n d e r . j p g ’ 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

) ; // SIVP t o o l b o x figure (1) ShowColorImage (I , ’ O r i g i n a l C o l o r Image ’ ) ; // IPD toolbox J = im2double ( I ) ; [ index , map ] = RGB2Ind ( I ) ; // IPD t o o l b o x pixels = prod ( size ( index ) ) ; hsv = rgb2hsv ( J ) ; h = hsv (: ,1) ; s = hsv (: ,2) ; v = hsv (: ,3) ; // F i n d s l o c a t i o n o f b l a c k and w h i t e p i x e l s darks = find (v <0.2) ; lights = find (s <0.05 & v >0.85) ; h ([ darks lights ]) = -1; // Gets t h e number o f a l l p i x e l s f o r e a c h c o l o u r b i n black_pixels = length ( darks ) / pixels ; white_pixels = length ( lights ) / pixels ; red = length ( find (( h > .9167 | h <= .083) & h ~= -1) ) / pixels ; yellow = length ( find ( h > .083 & h <= .25) ) / pixels ; 88

Figure 11.2: Read a Colour image and separate the colour image into red green and blue planes

89

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

green = length ( find ( h > .25 & h <= .4167) ) / pixels ; cyan = length ( find ( h > .4167 & h <= .5833) ) / pixels ; blue = length ( find ( h > .5833 & h <= .75) ) / pixels ; magenta = length ( find ( h > .75 & h <= .9167) ) / pixels ; // P l o t s h i s t o g r a m figure (2) a = gca () ; a . data_bounds =[0 ,0;8 ,1] n = 0:0.1:1; plot2d2 (n , red * ones (1 , length ( n ) ) ,5) n1 = 1:0.1:2; plot2d2 ( n1 , yellow * ones (1 , length ( n ) ) ,7) n2 = 2:0.1:3; plot2d2 ( n2 , green * ones (1 , length ( n ) ) ,8) n3 = 3:0.1:4; plot2d2 ( n3 , cyan * ones (1 , length ( n ) ) ,9) n4 = 4:0.1:5; plot2d2 ( n4 , blue * ones (1 , length ( n ) ) ,2) n5 = 5:0.1:6; plot2d2 ( n5 , magenta * ones (1 , length ( n ) ) ,3) n6 = 6:0.1:7; plot2d2 ( n6 , white_pixels * ones (1 , length ( n ) ) ,0) n7 = 7:0.1:8 plot2d2 ( n7 , black_pixels * ones (1 , length ( n ) ) ,5)

Scilab code Exa 11.18 Perform histogram equalisation of the given RGB image 1 2 3 4 5 6

// C a p t i o n : P e r f o r m h i s t o g r a m e q u a l i s a t i o n o f t h e g i v e n RGB image // F i g . 1 1 . 1 8 : MATLAB Example4 // p a g e 5 9 6 clc ; close ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 1 \ p e p p e r s . png ’ ) 90

7 8 9 10 11 12 13 14 15 16

; // SIVP t o o l b o x // c o n v e r s i o n o f RGB t o YIQ f o r m a t b = rgb2ntsc ( a ) ; // H i s t o g r a m e q u a l i s a t i o n o f Y component a l o n e b (: ,: ,1) = // c o n v e r s i o n o f YIQ t o RGB f o r m a t c = ntsc2rgb ( b ) ; figure (1) ShowColorImage (a , ’ O r i g i n a l Image ’ ) ; // IPD t o o l b o x figure (2) ShowColorImage (c , ’ H i s t o g t r a m e q u a l i z e d Image ’ ) ; // IPD t o o l b o x

Scilab code Exa 11.21 This program performs median filtering of the colour image 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

// C a p t i o n : T h i s program p e r f o r m s median f i l t e r i n g o f t h e c o l o u r image // F i g . 1 1 . 2 1 : MATLAB Example5 // p a g e 5 9 8 clc ; close ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 1 \ p e p p e r s . png ’ ) ; // SIVP t o o l b o x b = imnoise (a , ’ s a l t & p e p p e r ’ , 0.2) ; c (: ,: ,1) = MedianFilter ( b (: ,: ,1) , [3 3]) ; c (: ,: ,2) = MedianFilter ( b (: ,: ,2) , [3 3]) ; c (: ,: ,3) = MedianFilter ( b (: ,: ,3) , [3 3]) ; figure (1) ShowColorImage (a , ’ O r i g i n a l Image ’ ) ; // IPD t o o l b o x figure (2) ShowColorImage (b , ’ c o r r u p t e d Image ’ ) ; // IPD toolbox figure (3) ShowColorImage (c , ’ Median F i l t e r e d Image ’ ) ; // IPD 91

Figure 11.3: This program performs median filtering of the colour image toolbox

Scilab code Exa 11.24 Fitlering only the luminance component 1

// C a p t i o n : F i t l e r i n g o n l y t h e l u m i n a n c e component 92

2 // F i g . 1 1 . 2 4 : MATLAB Example6 3 // p a g e 5 9 9 4 clc ; 5 close ; 6 a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 1 \ p e p p e r s . png ’ ) 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

; // SIVP t o o l b o x // c o n v e r s i o n o f RGB t o YIQ f o r m a t yiq = rgb2ntsc ( a ) ; // E x t r a c t t h e Y component a l o n e b = yiq (: ,: ,1) ; h = [ -1 , -1 , -1; -1 ,8 , -1; -1 , -1 , -1]; // P e r f o r m h i g h p a s s f i l t e r i n g o n l y on Y component c1 = conv2d2 (b , h ) ; [m , n ]= size ( b ) ; for i =1: m for j =1: n D (i , j ) = c1 (i , j ) ; end end yiq (: ,: ,1) = D ; // c o n v e r t YIQ t o RGB f o r m a t a1 = ntsc2rgb ( yiq ) ; figure (1) ShowColorImage (a , ’ O r i g i n a l Image ’ ) ; // IPD t o o l b o x figure (2) ShowColorImage ( a1 , ’ High P a s s f i l t e r e d Image ’ ) ; // IPD t o o l b o x

Scilab code Exa 11.28 Perform gamma correction for the given colour image 1

// C a p t i o n : P e r f o r m gamma c o r r e c t i o n f o r t h e g i v e n c o l o u r image 93

Figure 11.4: Fitlering only the luminance component 2 // F i g . 1 1 . 2 8 : MATLAB Example7 3 // p a g e 6 0 3 4 close ; 5 clear all ; 6 clc ; 7 I = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 1 \ a r a r a u n a . png ’

);

// SIVP t o o l b o x

8 gamma_Value = 0.5; 9 max_intensity = 255; // f o r u i n t 8 image 10 // Look up t a b l e c r e a t i o n 11 LUT = max_intensity .*(([0: max_intensity ]./ 12 13 14 15 16 17 18 19 20 21

max_intensity ) .^ gamma_Value ) ; LUT = floor ( LUT ) ; // Mapping o f i n p u t p i x e l s i n t o l o o k u p t a b l e v a l u e s K = double ( I ) +1; J = zeros ( I ) ; [m ,n , p ]= size ( K ) ; for i = 1: m for j =1: n for k = 1: p J (i ,j , k ) = LUT ( K (i ,j , k ) ) ; end 94

Figure 11.5: Perform gamma correction for the given colour image 22 end 23 end 24 figure (1) 25 ShowColorImage (I , ’ O r i g i n a l Image ’ ) ; // IPD t o o l b o x 26 figure (2) 27 ShowColorImage ( uint8 ( J ) , ’Gamma C o r r e c t e d Image ’ ) ;

// IPD t o o l b o x

Scilab code Exa 11.30 Perform Pseudo Colouring Operation 1 // C a p t i o n : P e r f o r m Pseudo−C o l o u r i n g O p e r a t i o n 2 // F i g . 1 1 . 3 0 3 // p a g e 6 0 4 4 close ; 5 clear all ; 6 clc ; 7 K = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 1 \ l e n n a . j p g ’ ) ;

// SIVP t o o l b o x 8 [m , n ]= size ( K ) ; 9 I = uint8 ( K ) ;

95

10 for i = 1: m 11 for j =1: n 12 if ( I (i , j ) >=0 & I (i , j ) <50) 13 J (i ,j ,1) = I (i , j ) +50; 14 J (i ,j ,2) = I (i , j ) +100; 15 J (i ,j ,3) = I (i , j ) +10; 16 elseif ( I (i , j ) >=50 & I (i , j ) <100) 17 J (i ,j ,1) = I (i , j ) +35; 18 J (i ,j ,2) = I (i , j ) +128; 19 J (i ,j ,3) = I (i , j ) +10; 20 elseif ( I (i , j ) >=100 & I (i , j ) <150) 21 J (i ,j ,1) = I (i , j ) +152; 22 J (i ,j ,2) = I (i , j ) +130; 23 J (i ,j ,3) = I (i , j ) +15; 24 elseif ( I (i , j ) >=150 & I (i , j ) <200) 25 J (i ,j ,1) = I (i , j ) +50; 26 J (i ,j ,2) = I (i , j ) +140; 27 J (i ,j ,3) = I (i , j ) +25; 28 elseif ( I (i , j ) >=200 & I (i , j ) <=256) 29 J (i ,j ,1) = I (i , j ) +120; 30 J (i ,j ,2) = I (i , j ) +160; 31 J (i ,j ,3) = I (i , j ) +45; 32 end 33 end 34 end 35 figure (1) 36 ShowImage (K , ’ O r i g i n a l Image ’ ) ; // IPD t o o l b o x 37 figure (2) 38 ShowColorImage (J , ’ Pseudo C o l o u r e d Image ’ ) ; // IPD

toolbox

Scilab code Exa 11.32 Read an RGB image and segment it using the threshold method 96

Figure 11.6: Perform Pseudo Colouring Operation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

// C a p t i o n : Read an RGB image and s e g m e n t i t u s i n g t h e t h r e s h o l d method // F i g 1 1 . 3 2 // Page605 close ; clc ; I = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 1 \ a r a r a u n a . png ’ ) ; // SIVP t o o l b o x // C o n v e r s i o n o f RGB t o YCbCr b = rgb2ycbcr_1 ( I ) ; // SIVP t o o l b o x [m ,n , p ]= size ( b ) ; b = uint8 ( b ) ; // T h r e s h o l d i s a p p l i e d o n l y t o Cb component mask = b (: ,: ,2) >120; figure (1) ShowColorImage (I , ’ O r i g i n a l Image ’ ) ; // IPD t o o l b o x figure (2) ShowImage ( mask , ’ Segmented Image ’ ) ; // IPD t o o l b o x

97

Figure 11.7: Read an RGB image and segment it using the threshold method

98

Chapter 12 Wavelet based Image Processing

Scilab code Exa 12.9 Scilab code to perform wavelet decomposition 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

// C a p t i o n : S c i l a b c o d e t o p e r f o r m w a v e l e t decomposition // F i g 1 2 . 1 0 // Page624 clc ; close ; x = ReadImage ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 2 \ l e n n a . j p g ’ ); // The image i n u n s i g n e d i n t e g e r o r d o u b l e h a s t o be converted into normalized // d o u b l e f o r m a t x = im2double ( x ) ; // F i r s t L e v e l d e c o m p o s i t i o n [ CA , CH , CV , CD ]= dwt2 (x , ’ db1 ’ ) ; // S e c o n d l e v e l d e c o m p o s i t i o n [ CA1 , CH1 , CV1 , CD1 ]= dwt2 ( CA , ’ db1 ’ ) ; CA = im2int8 ( CA ) ; CH = im2int8 ( CH ) ; CV = im2int8 ( CV ) ; 99

17 18 19 20 21 22 23 24 25

CD = im2int8 ( CD ) ; CA1 = im2int8 ( CA1 ) ; CH1 = im2int8 ( CH1 ) ; CV1 = im2int8 ( CV1 ) ; CD1 = im2int8 ( CD1 ) ; A = [ CA , CH ; CV , CD ]; B = [ CA1 , CH1 ; CV1 , CD1 ]; imshow ( B ) title ( ’ R e s u l t o f S e c o n d L e v e l D e c o m p o s i t i o n ’ )

Scilab code Exa 12.42 Scilab code to generate different levels of a Gaussian pyramid 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

// C a p t i o n : S c i l a b c o d e t o g e n e r a t e d i f f e r e n t l e v e l s o f a G a u s s i a n pyramid // F i g 1 2 . 4 2 // Page651 clc ; close ; a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 2 \ a p p l e 3 . bmp ’ ) ; a = rgb2gray ( a ) ; b = a; kernelsize = input ( ’ E n t e r t h e s i z e o f t h e k e r n e l : ’ ) ; sd = input ( ’ E n t e r t h e s t a n d a r d d e v i a t i o n o f h t e G a u s s i a n window : ’ ) ; rf = input ( ’ E n t e r t h e R e d u c t i o n F a c t o r : ’ ) ; // R o u t i n e t o g e n e r a t e G a u s s i a n k e r n e l k = zeros ( kernelsize , kernelsize ) ; [ m n ] = size ( b ) ; t = 0; for i = 1: kernelsize for j =1: kernelsize k (i , j ) = exp ( -(( i - kernelsize /2) .^2+( j kernelsize /2) .^2) /(2* sd .^2) ) /(2* %pi * sd .^2) ; 100

19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

t = t + k (i , j ) ; end end for i = 1: kernelsize for j = 1: kernelsize k (i , j ) = k (i , j ) / t ; end end for t = 1:1: rf // c o n v o l v e i t w i t h t h e p i c t u r e FilteredImg = b ; if t ==1 FilteredImg = filter2 (k , b ) /255; else FilteredImg = filter2 (k , b ) ; end ; // compute t h e s i z e o f t h e r e d u c e d image m = m /2; n = n /2; // c r e a t e t h e r e d u c e d image t h r o u g h s a m p l i n g b = zeros (m , n ) ; for i = 1: m for j = 1: n b (i , j ) = FilteredImg ( i *2 , j *2) ; end ; end ; end ; figure ShowImage (a , ’ O r i g i n a l Image ’ ) figure ShowImage (b , ’ D i f f e r e n t L e v e l s o f G a u s a i n Pyramid ’ ) title ( ’ D i f f e r e n t L e v e l s o f G a u s a i n Pyramid L e v e l 2 ’ )

101

Figure 12.1: Scilab code to generate different levels of a Gaussian pyramid 102

Scilab code Exa 12.57 Scilab code to implement watermarking in spatial domain 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

// C a p t i o n : S c i l a b c o d e t o i m p l e m e n t w a t e r m a r k i n g i n s p a t i a l domain // F i g 1 2 . 5 7 // Page662 clc close a = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 2 \ cameraman . j p g ’ ); figure imshow ( a ) title ( ’ Base Image ’ ) ; b = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 2 \ k e y i m a g e . j p g ’ ); b = rgb2gray ( b ) ; b = imresize (b ,[32 32] , ’ b i c u b i c ’ ) ; [ m1 n1 ]= size ( b ) ; figure imshow ( b ) title ( ’ Mark Image ’ ) ; [ m n ]= size ( a ) ; i1 = 1; j1 = 1; p = 1; c = a; iii = 1; jjj = 1; a = uint8 ( a ) ; b = uint8 ( b ) ; for ff = 1:8 for i = 1:32 jjj = 1; for j = j1 : j1 + n1 -1 a (i , j ) = bitand ( a (i , j ) , uint8 (254) ) ; // LSB o f b a s e image i s s e t t o z e r o . temp = bitand ( b (i , jjj ) , uint8 ((2^ ff ) -1) ) ; 103

Figure 12.2: Scilab code to implement watermarking in spatial domain

32 33

//MSB o f t h e mark i s e x t r a c t e d . temp = temp /((2^ ff ) -1) ; c (i , j ) = bitor ( a (i , j ) , uint8 ( temp ) ) ; //MSB o f mark i s i n e r t e d i n t o t h e %LSB o f the base jjj = jjj +1;

34 35 end 36 end 37 j1 = j1 +32; 38 end 39 imshow ( c ) 40 title ( ’ Marked Image ’ ) ; 41 imwrite (c , ’E : \ DIP JAYARAMAN\ C h a p t e r 1 2 \ markimg . j p g ’ ) ;

Scilab code Exa 12.63 Scilab code to implement wavelet based watermarking 1

// C a p t i o n : S c i l a b c o d e t o i m p l e m e n t w a v e l e t −b a s e d watermarking 104

2 // F i g 1 2 . 6 3 3 // Page666 4 clc ; 5 close ; 6 // O r i g i n a l Image 7 img = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 2 \ cameraman .

jpg ’ ); 8 figure 9 imshow ( img ) 10 title ( ’ O r i g i n a l Image ’ ) ; 11 [ p q ] = size ( img ) ; 12 // G e n e r a t e t h e key 13 // key = i m r e a d ( ’ E : \ DIP JAYARAMAN\ C h a p t e r 1 2 \ keyimg1 .

png ’ ) ; 14 // key = i m r e s i z e ( key , [ p q ] ) ; 15 key = imread ( ’E : \ DIP JAYARAMAN\ C h a p t e r 1 2 \ k e y i m a g e . 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35

jpg ’ ); key = rgb2gray ( key ) ; c = 0.001; // I n i t i a l i s e t h e w e i g h t o f Watermarking figure imshow ( key ) title ( ’ Key ’ ) ; // Wavelet t r a n s f o r m o f o r i g i n a l image ( b a s e image ) img = double ( img ) ; key = double ( key ) ; [ ca , ch , cv , cd ] = dwt2 ( img , ’ db1 ’ ) ; // Compute 2D w a v e l e t transform // P e r f o r m t h e w a t e r m a r k i n g y = [ ca ch ; cv cd ]; Y = y + c * key ; p = p /2; q = q /2; for i =1: p for j =1: q nca (i , j ) = Y (i , j ) ; ncv (i , j ) = Y ( i +p , j ) ; nch (i , j ) = Y (i , j + q ) ; ncd (i , j ) = Y ( i +p , j + q ) ; 105

36 end 37 end 38 // D i s p l a y t h e Watermarked image 39 wimg = idwt2 ( nca , nch , ncv , ncd , ’ db1 ’ ) ; 40 wimg1 = uint8 ( wimg ) ; 41 figure 42 imshow ( wimg1 ) 43 title ( ’ Watermarked Image ’ ) 44 // E x t r a c t i o n o f key from Watermarked image 45 [ rca , rch , rcv , rcd ] = dwt2 ( wimg , ’ db1 ’ ) ; // Compute 2D 46 47 48 49 50 51 52

wavelet transform n1 =[ rca , rch ; rcv , rcd ]; N1 = n1 - y ; N1 = N1 *4; N1 = im2int8 ( N1 ) ; figure imshow ( N1 ) title ( ’ E x t r a c t t h e key from w a t e r m a r k e d image ’ )

106

Appendix Scilab code AP 1 2D Fast Fourier Transform 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

function [ a2 ] = fft2d ( a ) // a = any r e a l o r c o m p l e x 2D m a t r i x // a2 = 2D−DFT o f 2D m a t r i x ’a ’ m = size (a ,1) n = size (a ,2) // f o u r i e r t r a n s f o r m a l o n g t h e r o w s for i =1: n a1 (: , i ) = exp ( -2* %i * %pi *(0: m -1) ’.*.(0: m -1) / m ) * a (: , i ) end // f o u r i e r t r a n s f o r m a l o n g t h e c o l u m n s for j =1: m a2temp = exp ( -2* %i * %pi *(0: n -1) ’.*.(0: n -1) / n ) *( a1 (j ,:) ) .’ a2 (j ,:) = a2temp . ’ end for i = 1: m for j = 1: n if (( abs ( real ( a2 (i , j ) ) ) <0.0001) &( abs ( imag ( a2 ( i , j ) ) ) <0.0001) ) a2 (i , j ) =0; elseif ( abs ( real ( a2 (i , j ) ) ) <0.0001) a2 (i , j ) = 0+ %i * imag ( a2 (i , j ) ) ; elseif ( abs ( imag ( a2 (i , j ) ) ) <0.0001) a2 (i , j ) = real ( a2 (i , j ) ) +0; end end 107

25 end

Scilab code AP 2 2D Inverse FFT 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

function [ a ] = ifft2d ( a2 ) // a2 = 2D−DFT o f any r e a l o r c o m p l e x 2D m a t r i x // a = 2D−IDFT o f a2 m = size ( a2 ,1) n = size ( a2 ,2) // I n v e r s e F o u r i e r t r a n s f o r m a l o n g t h e r o w s for i =1: n a1 (: , i ) = exp (2* %i * %pi *(0: m -1) ’.*.(0: m -1) / m ) * a2 (: , i ) end // I n v e r s e f o u r i e r t r a n s f o r m a l o n g t h e c o l u m n s for j =1: m atemp = exp (2* %i * %pi *(0: n -1) ’.*.(0: n -1) / n ) *( a1 (j ,:) ) . ’ a (j ,:) = atemp . ’ end a = a /( m * n ) a = real ( a ) endfunction Scilab code AP 3 Median Filtering function

1 2 3 4 5 6 7 8 9 10 11 12 13

// The i n p u t t o t h e f u n c t i o n a r e t h e c o r r u p t e d image a and t h e d i m e n s i o n function [ Out_Imag ] = Func_medianall (a , N ) a = double ( a ) ; [ m n ]= size ( a ) ; Out_Imag = a ; if ( modulo (N ,2) ==1) Start =( N +1) /2; End = Start ; else Start = N /2; End = Start +1; end if ( modulo (N ,2) ==1) 108

14 limit1 =( N -1) /2; 15 limit2 = limit1 ; 16 else 17 limit1 =( N /2) -1; 18 limit2 = limit1 +1; 19 end 20 for i = Start :( m - End +1) , 21 for j = Start :( n - End +1) , 22 I =1; 23 for k = - limit1 : limit2 , 24 for l = - limit1 : limit2 , 25 mat ( I ) = a ( i +k , j + l ) ; 26 I = I +1; 27 end 28 end 29 mat = gsort ( mat ) ; // S o r t t h e e l e m e n t s t o 30 31 32 33 34 35 end 36 end

f i n d t h e median if ( modulo (N ,2) ==1) Out_Imag (i , j ) =( mat ((( N ^2) +1) /2) ) ; else Out_Imag (i , j ) =( mat (( N ^2) /2) + mat ((( N ^2) /2) +1) ) /2; end

Scilab code AP 4 To caculate gray level 1 function [ g ] = gray ( m ) 2 g = (0: m -1) ’/ max (m -1 ,1) 3 g = [g g g] 4 endfunction

Scilab code AP 5 To change the gray level of gray image 1 2

function [ bout ] = grayslice (I , z )

109

// Output v a r i a b l e s input v a r i a b l e s ) 4 bout =[]; 3

5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

i n i t i a l i s a t i o n ( not found in

// Number o f a r g u m e n t s i n f u n c t i o n c a l l [ %nargout , %nargin ] = argn (0) if %nargin ==1 then z = 10; elseif ~ type ( z ) ==1 then z = double ( z ) ; end ; n = z; if typeof ( I ) == ” u i n t 8 ” then z = (255*(0: n -1) ) / n ; elseif isa (I , ’ u i n t 1 6 ’ ) | isa (I , ’ i n t 1 6 ’ ) then z = 65535*(0:( n -1) ) / n ; else // I i s d o u b l e o r s i n g l e z = (0:( n -1) ) / n end ; [m , n ] = size ( I ) ; b = zeros (m , n ) ; // Loop o v e r a l l i n t e r v a l s , e x c e p t t h e l a s t for i = 1: length ( z ) -1 // j i s t h e i n d e x v a l u e we w i l l o u t p u t , s o i t depend upon s t o r a g e c l a s s if typeof ( b ) == ’ u i n t 8 ’ j = i -1; else j = i; end d = find (I >= z ( i ) & I < z ( i +1) ) ; if ~ isempty ( d ) , b(d) = j; end end // Take c a r e o f t h a t l a s t i n t e r v a l 110

39 d = find ( I >= z ( $ ) ) ; 40 if ~ isempty ( d ) then 41 // j i s t h e i n d e x v a l u e we w i l l o u t p u t , s o

depend upon s t o r a g e c l a s s 42 if typeof ( b ) == ” u i n t 8 ” then 43 j = length ( z ) -1; 44 else 45 j = length ( z ) ; 46 end ; 47 b(d) = j; 48 end ; 49 bout = b ; 50 bout = double ( bout ) ; 51 endfunction

111

it