Improving image clarity based on frequency filtering in Matlab

Introduction
To date, many algorithms have been developed to improve the quality of images characterized by the speed of the complexity of mathematical methods, the requirements for the resources of a computer system, etc. Moreover, 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 main idea of ​​filtering in the frequency domain is to select an accessory filter function that modifies F (u, v) in a specific way.
For example, a filter having an accessory function that, being multiplied by the centered function F (u, v), attenuates the high-frequency components of F (u, v), and leaves the low-frequency components relatively unchanged and is called the low-pass filter.
Images and their transformations are automatically considered periodic if filtering is implemented based on a two-dimensional Fourier transform.
When convolving periodic functions, interference effects arise between adjacent periods if these periods are located close to the length of the nonzero parts of the functions. Such interference is usually called a return error or an overlap error, which can be suppressed by supplementing the functions with zeros as follows:
Suppose that the functions f (x, y) and h (x, y) have dimensions AB and C D. We form two extended functions, both having dimensions 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)
  2. if nargin==1
  3.   PQ=2*AB;
  4. elseif nargin == 2&~ischar(CD)
  5.   PQ = AB+CD ;
  6.   PQ= 2 *ceil (PQ/2);
  7.   elseif nargin == 2
  8.     m= max (AB);
  9.     P=2^nextpow2(2*m);
  10.     PQ=[P, P];
  11. elseif nargin==3
  12.   m=max([AB CD]);
  13.    P=2^nextpow2(2*m);
  14.     PQ=[P, P];
  15. else
  16.   error('Неправильные входные данные');
  17. end
  18. 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 when the real and imaginary parts of the complex number are multiplied by the same number, the phase angle does not change. Such filters are commonly called filters with a zero phase shift.

  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


We generate a filter H in the frequency domain corresponding to the spatial Sobel 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 ('Low and high parameters 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 ('Invalid method.')
  33. end


The following dftuv function creates a grid array that is used in 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 that forms 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 is 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
  18.  


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 is 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




















Also popular now: