Image enhancement based on Matlab frequency filtering

Introduction
To date, many algorithms have been developed to improve the quality of images that differ in speed and complexity of mathematical methods, the requirements for computing system resources, etc. At the same time, one of the simplest methods is image processing based on its filtering in the frequency and spatial domains.

Basic concepts
The basis of linear filtering in the frequency and spatial domains is the convolution theorem. In fact, the basic idea of ​​filtering in the frequency domain is to select a filter accessory function that modifies F (u, v) in a specific way.
For example, a filter that has an ancillary function that, when multiplied by the centered function F (u, v), weakens the high-frequency components F (u, v) and leaves the low-frequency components relatively unchanged and is called a low-pass filter.
Images and their transformations are automatically considered periodic if the filtering is implemented on the basis of a two-dimensional Fourier transform.
When convoluting periodic functions, there are effects of interference between adjacent periods, if these periods are close relative to the length of the non-zero parts of the functions. Such an interference is usually called a return error or an overlap error, which can be suppressed by complementing the function with zeros as follows:
Suppose that the functions f (x, y) and h (x, y) have sizes AB and C D. We form two extended functions, both having sizes P and Q, by adding zeros to f and g.
P≥A + C + 1 and Q≥B + D + 1
We calculate the minimum even values ​​of P and Q that satisfy these conditions.
We implement the above algorithm as a function.

1. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
2. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
3. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
4. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
5. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
6. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
7. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
8. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
9. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
10. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
11. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
12. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
13. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
14. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
15. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
16. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
17. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .
18. function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .

function PQ = paddedsize(AB, CD, PARAM) if nargin==1 PQ=2*AB; elseif nargin == 2&~ischar(CD) PQ = AB+CD ; PQ= 2 *ceil (PQ/2); elseif nargin == 2 m= max (AB); P=2^nextpow2(2*m); PQ=[P, P]; elseif nargin==3 m=max([AB CD]); P=2^nextpow2(2*m); PQ=[P, P]; else error ( ' '); end end * This source code was highlighted with Source Code Highlighter .

The transfer function of the filter H (u, v) is multiplied by the real and imaginary parts of F (u, v). If the function H (u, v) was real, then the phase part of the product does not change, as can be seen from the phase equation, since multiplying the real and imaginary parts of the complex number by the same number does not change the phase angle. Such filters are called zero phase shift filters.
')

1. function g = dftfilt ( f, H )
2. F = fft2 ( f, size ( H, 1 ) , size ( H, 2 ) ) ;
3. Gi = H. * F;
4. g = real ( ifft2 ( Gi ) ) ;
5. g = g ( 1 : size ( f, 1 ) , 1 : size ( f, 2 ) ) ;
6. end

Generate the H filter in the frequency domain corresponding to the Sobel spatial filter, which improves the vertical edges of the image.

1. function g = gscale ( f, varargin )
2. if length ( varargin ) == 0
3. method = 'full8' ;
4. else
5. method = varargin { 1 } ;
6. end
7. if strcmp ( class ( f ) , 'double' ) & ( max ( f ( : ) ) > 1 | min ( f ( :) ) < 0 )
8. f = mat2gray ( f ) ;
9. end
10. switch method
11. case 'full8'
12. g = im2uint8 ( mat2gray ( double ( f ) ) ) ;
13. case 'full16'
14. g = im2uint16 ( mat2gray ( double ( f ) ) ) ;
15. case 'minmax'
16. low = varargin { 2 } ; high = varargin { 3 } ;
17. if low> 1 | low < 0 | high> 1 | high < 0
18. error ( 'Parameters low and high must be changed' )
19. end
20. if strcmp ( class ( f ) , 'double' )
21. low_in = min ( f ( :) ) ;
22. high_in = max ( f ( :) ) ;
23. elseif strcmp ( class ( f ) , 'uint8' )
24. low_in = min ( f ( :) ) ./ 255 ;
25. high_in = max ( f ( :) ) ./ 255 ;
26. elseif strcmp ( class ( f ) , 'uint16' )
27. low_in = min ( f ( :) ) ./ 65535 ;
28. high_in = max ( f ( :) ) ./ 65535 ;
29. end
30. g = imadjust ( f, [ low_in, high_in ] , [ low, high ] ) ;
31. otherwise
32. error ( 'Wrong method.' )
33. end

The following dftuv function creates a grid array, which is used when calculating distances and other similar actions.

1. function [ U, V ] = dftuv ( M, N )
2. u = 0 : ( M ) ;
3. v = 0 : ( N ) ;
4. idx = find ( u> M / 2 ) ;
5. u ( idx ) = u ( idx ) ;
6. idy = find ( v> N / 2 ) ;
7. v ( idy ) = v ( idy ) ;
8. [ V, U ] = meshgrid ( v, u ) ;
9. end

I will give an example of creating a function forming a low-pass filter.

1. function [ H, D ] = lpfilter ( type , M, N, D0, n )
2. [ U, V ] = dftuv ( M, N ) ;
3. D = sqrt ( U. ^ 2 + V. ^ 2 ) ;
4. switch type
5. case 'ideal'
6. H = double ( D <= D0 ) ;
7. case 'btw'
8. if nargin == 4
9. n = 1 ;
10. end
11. H = 1. / ( 1 + ( D / D0 ) . ^ ( 2 * n ) ) ;
12. case 'gaussian'
13. H = exp ( ( D. ^ 2 ) ./ ( 2 * ( D0 ^ 2 ) ) ) ;
14. otherwise
15. error ( 'Unknown filter type' ) ;
16. end
17. end

Now an example of creating a function forming a high-pass filter.

1. function H = hpfilter ( type , M, N, D0, n )
2. [ U, V ] = dftuv ( M, N ) ;
3. D = sqrt ( U. ^ 2 + V. ^ 2 ) ;
4. switch type
5. case 'ideal'
6. H = double ( D <= D0 ) ;
7. case 'btw'
8. if nargin == 4
9. n = 1 ;
10. end
11. H = 1 - ( 1. / ( 1 + ( D / D0 ) . ^ ( 2 * n ) ) ) ;
12. case 'gaussian'
13. H = 1 - ( exp ( ( D. ^ 2 ) ./ ( 2 * ( D0 ^ 2 ) ) ) ) ;
14. otherwise
15. error ( 'Unknown filter type' ) ;
16. end
17. end

And the script itself to increase the clarity of images, using the above functions.

1. img = imread ( 'D: \ Photo \ nature.jpg' ) ;
2. red = img ( :,:, 1 ) ;
3. F = fft2 ( img ) ;
4. S = fftshift ( log ( 1 + abs ( F ) ) ) ;
5. S = gscale ( S ) ;
6. imshow ( img ) , figure , imshow ( S ) ;
7. PQ = paddedsize ( size ( red ) ) ;
8. [ U, V ] = dftuv ( PQ ( 1 ) , PQ ( 2 ) ) ;
9. D0 = 0.05 * PQ ( 2 ) ;
10. F = fft2 ( red, PQ ( 1 ) , PQ ( 2 ) ) ;
11. H = exp ( - ( U. ^ 2 + V. ^ 2 ) / ( 2 * ( D0 ^ 2 ) ) ) ;
12. g = dftfilt ( red, H ) ;
13. figure , imshow ( fftshift ( H ) , [ ] ) ;
14. figure , mesh ( H ( 1 : 10 : 500 , 1 : 10 : 500 ) )
15. axis ( [ 0 50 0 50 0 1 ] )
16. colormap ( gray )
17. grid off
18. axis off
19. view ( - 25 , 0 )
20. H = fftshift ( hpfilter ( 'gaussian' , 500 , 500 , 50 ) ) ;
21. mesh ( H ( 1 : 10 : 500 , 1 : 10 : 500 ) )
22. axis ( [ 0 50 0 50 0 1 ] )
23. colormap ( [ 0 0 0 ] )
24. grid off
25. axis off
26. view ( - 163 , 64 )
27. figure , imshow ( H, [ ] ) ;
28. PQ = paddedsize ( size ( red ) ) ;
29. D0 = 0.05. * PQ ( 1 ) ;
30. H = hpfilter ( 'gaussian' , PQ ( 1 ) , PQ ( 2 ) , D0 ) ;
31. g = dftfilt ( red, H ) ;
32. figure , imshow ( g, [ 0 255 ] ) ;
33. PQ = paddedsize ( size ( red ) ) ;
34. D0 = 0.05. * PQ ( 1 ) ;
35. HBW = hpfilter ( 'btw' , PQ ( 1 ) , PQ ( 2 ) , D0, 2 ) ;
36. H = 0.5 + 2 * HBW;
37. gbw = dftfilt ( red, HBW ) ;
38. gbw = gscale ( gbw ) ;
39. ghf = dftfilt ( red, H ) ;
40. gbf = gscale ( ghf ) ;
41. img1 = histeq ( red, 256 ) ;
42. ghe = histeq ( ghf, 256 ) ;
43. figure , imshow ( red ) ;
44. figure , imshow ( ghe, [ ] ) ;
45. figure , imshow ( img1, [ ] ) ;

**** Examples of work ****


















Final result