Фильтр High Pass Butterworth на изображениях в MATLAB

Question

Мне нужно реализовать фильтр Баттерворта с высоким разрешением в MATLAB для целей фильтрации изображений. Я реализовал один, но похоже, что он не работает. Вот код, который я написал. Может ли кто-нибудь сказать мне, что не так?

n=1; 
d=50; 
A=1.5; 
im=imread('imagex.jpg'); 
h=size(im,1); 
w=size(im,2); 
[x y]=meshgrid(-floor(w/2):floor(w-1/2),-floor(h/2):floor(h-1/2)); 
hhp=(1./(d./(x.^2+y.^2).^0.5).^(2*n)); 
image_2Dfilter=fftshift(fft2(im)); 
Image_butterworth=image_2Dfilter; 
imshow(Image_butterworth); 
ifftshow(Image_butterworth); 

Answer 1

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

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

Вы можете найти величину спектра с помощью функции abs. Даже если вы это сделаете, если вы сделали imshow непосредственно на величину, вы получите визуализацию, равную нулю всюду за исключением для середины. Это связано с тем, что компонент постоянного тока настолько велик, а остальная часть спектра мала для сравнения.

Позвольте мне показать вам пример. Это кинооператор изображение, которое является частью обработки изображений инструментов:

im = imread('cameraman.tif'); 
figure; 
imshow(im); 

Теперь давайте визуализировать спектр и обеспечение того, чтобы компонент постоянного тока находится в центре изображения - вы уже сделали это с fftshift. Также неплохо, чтобы лидировал изображение до double, чтобы обеспечить максимальную точность данных. Кроме того, убедитесь, что вы применяете abs найти величину:

fftim = fftshift(fft2(double(im))); 
mag = abs(fftim); 
figure; 
imshow(mag, []); 

Как вы можете видеть, это не очень полезно по той причине, о которой я говорил. Лучшим способом визуализации спектра изображения обычно является применение преобразования log в спектр. Это также полезно, если вы хотите уменьшить или удалить среднее значение, чтобы динамический диапазон лучше подходит для отображения. Другими словами, вы должны добавить 1 к величине, а затем применить логарифм к величине, чтобы более высокие значения могли уменьшаться. Это не имеет значения, какую базу вы используете, так что я буду просто использовать натуральный логарифм, который инкапсулируется с помощью команды log:

figure; 
imshow(log(1 + mag), []); 

Теперь это много лучше. Теперь мы перейдем к вашему механизму фильтрации. Ваш фильтр Баттерворта немного некорректен. Координаты meshgrid немного ошибочны.Операция -1 это на окончание интервала должен выходить за пределы:

[x y]=meshgrid(-floor(w/2):floor(w/2)-1,-floor(h/2):floor(h/2)-1); 

Помните, что вы определяете симметричный интервал относительно центра изображения, и то, что у вас изначально не было правильным. Я также хотел бы упомянуть, что это выглядит как фильтр верхних частот, поэтому выход должен выглядеть как обнаружение краев. Кроме того, неверно определение фильтра высоких частот Баттерворта. Правильное определение фильтра в частотной области:

D(u,v) является расстоянием от центра изображения в частотной области, Do является расстоянием отсечки, а B является контролирующим масштабным коэффициентом контролируя то, что желаемая коэффициент усиления будет на расстоянии отсечки. n - это порядок фильтра. Do в вашем случае это d = 50. На практике B = sqrt(2) - 1 так, что на расстоянии отсечки Do, D(u,v) = 1/sqrt(2) = 0.707, которое является частотой отсечки 3 дБ, чаще всего наблюдаемой в электронных фильтрах цепи. Иногда вы можете увидеть, что B установлено в 1 для простоты, но обычно устанавливается значение B = sqrt(2) - 1.

Однако ваш текущий код не выполняет фильтрацию. Чтобы фильтровать в частотной области, вы просто умножаете спектр изображения на спектр самого фильтра. Это эквивалентно свертке в пространственной области. Как только вы это сделаете, вы просто отмените fftshift, который был выполнен на изображении, возьмите обратный БПФ и затем устраните любые мнимые компоненты, которые вызваны численной неточностями. Кроме того, давайте перейдем к uint8, чтобы убедиться, что мы уважаем исходный тип изображения.

Это можно сделать так:

%// Your code with meshgrid fix 
n=1; 
d=50; 
h=size(im,1); 
w=size(im,2); 
fftim = fftshift(fft2(double(im))); 
[x y]=meshgrid(-floor(w/2):floor(w/2)-1,-floor(h/2):floor(h/2)-1); 
%hhp=(1./(d./(x.^2+y.^2).^0.5).^(2*n)); 

%%%%%%// New code 
B = sqrt(2) - 1; %// Define B 
D = sqrt(x.^2 + y.^2); %// Define distance to centre 
hhp = 1 ./ (1 + B * ((d ./ D).^(2 * n))); 
out_spec_centre = fftim .* hhp; 

%// Uncentre spectrum 
out_spec = ifftshift(out_spec_centre); 

%// Inverse FFT, get real components, and cast 
out = uint8(real(ifft2(out_spec))); 

%// Show image 
imshow(out); 

Если вы хотите увидеть, что отфильтрованный спектр выглядит, просто сделать это:

figure; 
imshow(log(1 + abs(out_spec_centre)), []); 

Получаем:

Это имеет смысл. Вы видите, что в середине спектра он немного темнее по сравнению с внешними краями спектра. Это потому, что с фильтром Баттерворта с высокими пропусками вы усиливаете более высокие частотные термины и визуализируется как более высокая интенсивность.

Теперь out содержит ваше обработанное изображение, и мы, наконец, получить это:

Это выглядит как прекрасный результат! Однако наивное отбрасывание изображения до uint8 усекает любые отрицательные значения до 0 и любые положительные значения, превышающие 255-255. Поскольку это обнаружение краев, вы хотите обнаружить как отрицательный, так и положительный переходы ... так что хорошая идея была бы до нормализовать выход так, чтобы он варьировался от [0,1], а затем с uint8 после умножения на 255. Таким образом, никакие изменения в изображении не будут визуализироваться серым цветом, негативные изменения визуализируются как темные, а положительные изменения визуализируются как белые ....так что вы хотите сделать что-то вроде этого:

%// Your code with meshgrid fix 
n=1; 
d=50; 
h=size(im,1); 
w=size(im,2); 
fftim = fftshift(fft2(double(im))); 
[x y]=meshgrid(-floor(w/2):floor(w/2)-1,-floor(h/2):floor(h/2)-1); 
%hhp=(1./(d./(x.^2+y.^2).^0.5).^(2*n)); 

%%%%%%// New code 
B = sqrt(2) - 1; %// Define B 
D = sqrt(x.^2 + y.^2); %// Define distance to centre 
hhp = 1 ./ (1 + B * ((d ./ D).^(2 * n))); 
out_spec_centre = fftim .* hhp; 

%// Uncentre spectrum 
out_spec = ifftshift(out_spec_centre); 

%// Inverse FFT, get real components 
out = real(ifft2(out_spec)); 

%// Normalize and cast 
out = (out - min(out(:)))/(max(out(:)) - min(out(:))); 
out = uint8(255*out); 

%// Show image 
imshow(out); 

Мы получаем это:

Answer 2

Я думаю, что вы должны работать немножко Diferent

n=1; 
D0=50; % change the name for d0, d is usuaally the (u²+v²)⁽1/2) 

A=1.5; % normally the amplitude is 1 

im=imread('cameraman.jpg'); 

[M,N]=size(im); % is easy to get the h and w like this 

% compute the 2d fourier transform in order to multiply 

F=fft2(double(im)); 

% compute your filter and do the meshgrid for your matrix but it is M*n, and get only the real part 

u=0:(M-1); 
v=0:(N-1); 

idx=find(u>M/2); 
u(idx)=u(idx)-M; 
idy=find(v>N/2); 
v(idy)=v(idy)-N; 
[V,U]=meshgrid(v,u); 
D=sqrt(U.^2+V.^2); 
H =A * (1./(1 + (D0./D).^(2*n))); 

% multiply element by element 

G=H.*F; 
g=real(ifft2(double(G))); 
subplot(1,2,1); imshow(im); title('Input image'); 
subplot(1,2,2); imshow(g,[ ]); title('filtered image');