Как стать автором
Обновить

Повышение четкости изображений на основе частотной фильтрации в Matlab

Обработка изображений *Matlab *
Из песочницы
Введение
На сегодняшний день разработано множество алгоритмов по улучшению качества изображений отличающихся быстродействием сложностью математических методов требованиям к ресурсам вычислительной системы и т.д. При этом одним из наиболее простых методов является обработка изображений на основе его фильтрации в частотной и пространственной областях.

Базовые концепции
Основой линейной фильтрации в частотной и пространственной области служит теорема о свертке. На самом деле, основная идея фильтрации в частотной области заключается в подборе придаточной функции фильтра, которая модифицирует F(u ,v), специфическим образом.
Например, фильтр, имеющий придаточную функцию, которая будучи умноженной на центрированную функцию F(u ,v), ослабляет высокочастотные компоненты F(u ,v), и оставляет низкочастотные компоненты относительно неизменными и называется низкочастотным фильтром.
Изображения и их преобразования автоматически считаются периодическими, если фильтрация реализована на основе двумерного преобразования Фурье.
При свертки периодических функций, возникают эффекты интерференции между смежными периодами, если эти периоды расположены близко относительно длины ненулевых частей функций. Такую интерференцию принято называть ошибкой возврата или ошибкой перекрытия, которую можно подавить дополняя функции нулями следующим образом:
Предположим что функции f(x, y) и h(x, y) имеют размеры A B и C D. Сформируем две расширенные функции, обе имеющие размеры P и Q, путем добавления нулей к f и g.
P≥A+C+1 и Q≥B+D+1
Вычислим минимальные четные значения P и Q которые удовлетворяют этим условиям.
Реализуем вышеописанный алгоритм как функцию.

  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
.


Передаточная функция фильтра H(u, v) умножается на вещественную и мнимую части F(u, v). Если функция H(u, v) была, вещественной то фазовая часть произведения не меняется, что видно из фазового уравнения, так как при умножении вещественной и мнимой части комплексного числа на одно и то же число фазовый угол не меняется. Такие фильтры принято называть фильтрами с нулевым сдвигом фазы.

  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


Сгенерируем фильтр H в частотной области, соответствующий пространственному фильтру Собела, который улучшает вертикальные края изображения.

  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(:))>1min (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 и high должны быть изменены')
  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 ('Неправильный метод.')
  33. end


Следущая функция dftuv создает сеточный массив, который используется при вычислении расстояний и при других подобных действиях.

  1. function [U , V]= dftuv( M, N )
  2. u=0:();
  3. v=0:();
  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


Приведу пример создания функции формирующей низкочастотный фильтр.

  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 ('Неизвестный тип фильтра');
  16.     end
  17. end
  18.  


Теперь пример создания функции формирующей высокочастотный фильтр.

  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 ('Неизвестный тип фильтра');
  16.     end
  17. end


И сам скрипт для увеличения четкости изображений, с использованием вышеописанных функций.

  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. figuremesh(H(1:10:5001:10:500))
  15. axis([0 50 0 50 0 1])
  16. colormap(gray)
  17. grid off
  18. axis off
  19. view (-250)
  20. H=fftshift (hpfilter ('gaussian'50050050));
  21. mesh(H(1:10:5001:10:500))
  22. axis([0 50 0 50 0 1])
  23. colormap([0 0 0])
  24. grid off
  25. axis off
  26. view (-16364)
  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, []);


****Примеры работы****


















Конечный результат
Теги: четкость изображенийчастотная фильтрацияmatlab
Хабы: Обработка изображений Matlab
Всего голосов 7: ↑5 и ↓2 +3
Комментарии 8
Комментарии Комментарии 8

Похожие публикации

Лучшие публикации за сутки