Создание фильтра высоких частот в matlab

Question

Создание фильтра высоких частот в matlab

Я пытаюсь создать фильтр высоких частот в Matlab. Я генерирую гауссово ядро, используя

function kernel = compute_kernel(sigma,size)

[x,y] = meshgrid(-size/2:size/2,-size/2:size/2);

constant = 1/(2*pi*sigma*sigma);

kernel = constant*exp( -(y.^2 + x.^2 )/(2 * sigma * sigma));

kernel = (kernel - min(kernel(:)))./(max(kernel(:)) - min(kernel(:)));

end

Затем после создания ядра я использую его для создания фильтра низких частот для изображения (переменная im2 ):

g = compute_kernel(9,101);

im2_low = conv2(im2,g,'same');

Как я понимаю, я могу затем использовать вычитание отфильтрованного изображения из исходного изображения (в частотной области) для извлечения высоких частот, что делает его эквивалентом фильтра высоких частот.

F = fft2(im2_low);

IM2 = fft2(im2);

IM2_high = IM2 - F;



figure; fftshow(IM2_high);

im2_high = ifft2(IM2_high);

figure; imshow(im2_high,[]);

Кажется, что-то здесь не так. Когда я просмотр высокочастотного отфильтрованного изображения он кажется цветным инвертированным размытым изображением, а не тем, с краями, определенными, как я видел в интернете. Я не уверен, является ли мой процесс неправильным или я просто использую неправильные значения для моего ядра Гаусса.

902 2

image-processing matlab

2 ответов:

Comments

Ничего не найдено.

Ander Biguri · Accepted Answer · 2018-01-16 10:56:48

Любое ядро, которое вы хотите использовать для поддержания функций изображения (то есть вы не хотите, чтобы величина чего-то, но изображение выглядело как узнаваемый образ человека), вы должны убедиться, что вы делаете что-то с ядром: нормализуете его.
Вы, кажется, пытались это сделать, но неправильно истолковали значение нормализации в ядрах. Не нужно быть [0-1], их сумма должна быть 1.
Итак, берем ваш код:
im2=imread('https://upload.wikimedia.org/wikipedia/en/2/24/Lenna.png');
im2=double(rgb2gray(im2));

sigma=9;
sizei=101;
[x,y] = meshgrid(-sizei/2:sizei/2,-sizei/2:sizei/2);
constant = 1/(2*pi*sigma*sigma);
kernel = constant*exp( -(y.^2 + x.^2 )/(2 * sigma * sigma));
%%%%%% NORMALIZATION
kernel=kernel/sum(kernel(:));
%%%%%%

im2_low = conv2(im2,kernel,'same');

F = fft2(im2_low);
IM2 = fft2(im2);
IM2_high = IM2 - F;


im2_high = ifft2(IM2_high);
figure; imshow(im2_high,[]);
Но, как Крислуэнго упоминает, что вычитание-это операция, которая не изменяется в области Фурье, поэтому ответ
im2_high=im2-im2_low

Cris Luengo · Accepted Answer · 2018-01-17 04:08:55

Это длинный ответ на короткий вопрос. Прочтите его, если хотите что-то узнать.
Фильтр нижних частот и фильтр верхних частот являются линейными фильтрами. Линейный фильтр может быть применен в пространственной области через свертку или в частотной области (она же область Фурье) в качестве умножения.
Верно, что в области Фурье различие между ядром фильтра нижних частот и тождественным фильтром (all-pass filter) является высокочастотным фильтром. фильтр:
high_pass_filter = identity_filter - low_pass_filter
Фильтр идентификации будет ядром, где каждый элемент равен 1. Фильтр применяется путем умножения, поэтому
IM2 * high_pass_filter = IM2 * ( identity_filter - low_pass_filter )
Что то же самое, что
IM2 * high_pass_filter = IM2 - IM2 * low_pass_filter
(здесь, как и в вопросе, IM2 является представлением изображения в области Фурье im2; все элементы с желтыми границами предназначены для уравнений, но записываются в псевдокоде, с символом *, используемым для умножения).

Таким образом, ОП хочет применить фильтр низких частот и вычесть входное изображение в области Фурье, чтобы получить изображение с фильтром высоких частот.
Однако одно из свойств преобразования Фурье состоит в том, что оно является линейным преобразованием. Это означает, что

F(a*f + b*g) == a * F(f) + b * F(g)

(с F(.) преобразованием Фурье, a и b константами, а также f и g функциями). Устанавливая a=1 и b=-1, а также g низкочастотное фильтрованное изображение и f входное изображение, мы get

F(im2 - im2_low) == F(im2) - F(im2_low)

То есть вычитание в пространственной области и в области Фурье эквивалентны. Таким образом, если вычислить im2_low в пространственной области, то нет необходимости переходить в область Фурье для вычитания. Эти два бита кода дают идентичный результат (с точностью до цифр):
F = fft2(im2_low); IM2 = fft2(im2); IM2_high = IM2 - F; im2_high = ifft2(IM2_high); im2_high = im2 - im2_low;
Кроме того, свертка также линейна. Это означает, что, если вы думаете о F(.) в уравнениях выше как о свертке, эти уравнения все еще остаются в силе. Вы может делать такие манипуляции:

conv(f, h) - f == conv(f, h) - conv(f, 1) == conv(f, h-1)

Это непосредственно приводит к определению фильтра высоких частот в пространственной области:

g = - compute_kernel(9,101); g(51,51) = g(51,51) + 1; im2_high2 = conv2(im2,g,'same');

Вы увидите, что max(max(abs(im2_high-im2_high2))) дает значение, очень близкое к 0.

Замечание относительно вычисления фильтра Гаусса:

Функция compute_kernel, опубликованная в вопросе, вычисляет ядро 2D-фильтра, непосредственно оценивая 2D-Гаусса. Результирующее ядро фильтра составляет 101x101 пикселов, что означает, что вычисление свертки требует 101 * 101 * N умножения и сложения (MADs) с N числом пикселей в отфильтрованном изображении. Однако гауссов фильтр является сепарабельным, что означает, что тот же результат может быть получен только в 101 * 2 * N Мадах (в 50 раз меньше!). Кроме того, для sigma = 9 можно обойтись и меньшим ядром.

Размер ядра Гаусса:
Гауссова функция никогда не достигает нуля, но она очень быстро приближается к нулю. Когда режешь его в 3*sigma, очень мало из этого теряется. Я считаю, что 3 Сигмы-это хороший баланс. В случае sigma = 9 отсечение 3 Сигмы приводит к ядру с 55 пикселями (3*sigma * 2 + 1).

Гауссова отделимость:

Многомерного Гаусса можно получить, умножив 1D Гауссианов вместе:

exp(-(y.^2+x.^2)/(2*sigma*sigma)) == exp(-(x.^2)/(2*sigma*sigma)) * exp(-(y.^2)/(2*sigma*sigma))

Это приводит к гораздо более эффективной реализации свертки:

conv(f,h1*h2) == conv( conv(f,h1), h2 )

То есть свертывание образа с помощью фильтр столбцов h1, а затем свертка результата с помощью фильтра строк h2 аналогична свертке изображения с помощью 2D-фильтра h1*h2. В коде:

sigma = 9; sizei = ceil(3*sigma); % 3 sigma cutoff g = exp(-(-sizei:sizei).^2/(2*sigma.^2)); % 1D Gaussian kernel g = g/sum(g(:)); % normalize kernel im2_low = conv2(g,g,im2,'same'); g2d = g' * g; im2_low2 = conv2(im2,g2d,'same');

Разница в числовой неточности:

max(max(abs(im2_low-im2_low2)))

ans = 1.3927e-12

Вы найдете более подробное описание гауссовой фильтрации в моем блоге, а также некоторые проблемы, с которыми вы можете столкнуться при использовании Matlab's Image Processing Toolbox.