GoooodMan 16 сен 2011 в 13:39

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

4 мин

16K

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

Базовые концепции
Основой линейной фильтрации в частотной и пространственной области служит теорема о свертке. На самом деле, основная идея фильтрации в частотной области заключается в подборе придаточной функции фильтра, которая модифицирует 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 которые удовлетворяют этим условиям.
Реализуем вышеописанный алгоритм как функцию.

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.

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

function g =dftfilt(f, H)

F=fft2(f, size(H,1), size(H,2));

Gi=H.*F;

g=real(ifft2(Gi));

g=g(1:size(f,1),1:size(f,2));

end

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

function g = gscale( f, varargin )

if length(varargin)==0

    method='full8';

else

    method = varargin{1};

end

if strcmp(class(f), 'double')&(max(f(:))>1| min (f(:))<0)

    f=mat2gray(f);

end

switch method

    case 'full8'

       g=im2uint8(mat2gray(double (f)));

    case 'full16'

       g=im2uint16(mat2gray(double (f)));

    case 'minmax'

       low=varargin{2}; high=varargin{3};

       if low>1 | low<0 | high>1 | high<0

           error('Параметры low и high должны быть изменены')

       end

       if strcmp(class(f), 'double')

           low_in=min(f(:));

           high_in=max(f(:));

       elseif strcmp(class(f), 'uint8')

           low_in=min(f(:))./255;

           high_in=max(f(:))./255;

       elseif strcmp(class(f), 'uint16')

           low_in=min(f(:))./65535;

           high_in=max(f(:))./65535;

       end

       g=imadjust(f, [low_in, high_in], [low, high]);

    otherwise

        error ('Неправильный метод.')

end

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

function [U , V]= dftuv( M, N )

u=0:(M );

v=0:(N );

idx =find(u>M/2);

u(idx) = u(idx);

idy=find(v>N/2);

v(idy) = v(idy) ;

[V, U]=meshgrid(v, u);

end

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

function [H, D] = lpfilter( type, M, N, D0, n )

[U, V]=dftuv(M, N);

D=sqrt(U.^2+V.^2);

switch type

    case 'ideal'

        H=double (D<=D0);

    case 'btw'

        if nargin==4

            n=1;

        end

        H=1./(1+(D/D0).^(2*n));

    case 'gaussian'

        H= exp ((D.^2)./(2*(D0^2)));

    otherwise

        error ('Неизвестный тип фильтра');

    end

end

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

function H = hpfilter( type, M, N, D0, n )

   [U, V]=dftuv(M, N);

D=sqrt(U.^2+V.^2);

switch type

    case 'ideal'

        H=double (D<=D0);

    case 'btw'

        if nargin==4

            n=1;

        end

        H=1-(1./(1+(D/D0).^(2*n)));

    case 'gaussian'

        H=1- (exp ((D.^2)./(2*(D0^2))));

    otherwise

        error ('Неизвестный тип фильтра');

    end

end

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

img = imread('D:\Photo\nature.jpg');

red = img(:,:, 1);

F=fft2(img);

S=fftshift(log(1+abs(F)));

S=gscale(S);

imshow(img), figure, imshow (S);

PQ=paddedsize(size(red));

[U, V]=dftuv(PQ(1), PQ(2));

D0=0.05*PQ(2);

F=fft2(red, PQ(1), PQ(2));

H=exp(-(U.^2+V.^2)/(2*(D0^2)));

g=dftfilt(red, H);

figure, imshow(fftshift(H), []);

figure, mesh(H(1:10:500, 1:10:500))

axis([0 50 0 50 0 1])

colormap(gray)

grid off

axis off

view (-25, 0)

H=fftshift (hpfilter ('gaussian', 500, 500, 50));

mesh(H(1:10:500, 1:10:500))

axis([0 50 0 50 0 1])

colormap([0 0 0])

grid off

axis off

view (-163, 64)

figure, imshow (H, []);

PQ=paddedsize (size(red));

D0=0.05.*PQ(1);

H=hpfilter('gaussian', PQ(1), PQ(2), D0);

g=dftfilt(red, H);

figure, imshow (g, [0 255]);

PQ=paddedsize (size(red));

D0=0.05.*PQ(1);

HBW=hpfilter('btw', PQ(1), PQ(2), D0,2);

H=0.5+2*HBW;

gbw=dftfilt(red, HBW);

gbw=gscale(gbw);

ghf=dftfilt(red, H);

gbf=gscale(ghf);

img1=histeq(red, 256);

ghe=histeq(ghf, 256);

figure, imshow (red);

figure, imshow (ghe, []);

figure, imshow (img1, []);