Как применить всплески Габора к изображению?

Ответ 1

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

Вы можете сверлить фильтр с изображением, используя теорему свертки, взяв обратное преобразование Фурье элементарного произведения преобразований Фурье изображения и фильтр. Вот основная формула:

%# Our image needs to be 2D (grayscale)
if ndims(img) > 2;
    img = rgb2gray(img);
end
%# It is also best if the image has double precision
img = im2double(img);

[m,n] = size(img);
[mf,nf] = size(GW);
GW = padarray(GW,[n-nf m-mf]/2);
GW = ifftshift(GW);
imgf = ifft2( fft2(img) .* fft2(img) );

Как правило, свертка FFT превосходит ядра размером > 20. Для получения дополнительной информации я рекомендую Numerical Recipes в C, который имеет хорошее языковое агностическое описание метода и его оговорки.

Ваши ядра уже большие, но с помощью метода FFT они могут быть такими же большими, как и изображение, поскольку они заполняются до такого размера независимо. Из-за периодичности БПФ метод выполняет круговую свертку. Это означает, что фильтр будет обтекать границы изображения, поэтому нам нужно также проложить изображение, чтобы устранить этот краевой эффект. Наконец, поскольку мы хотим, чтобы общий отклик всем фильтрам (по крайней мере, в типичной реализации), мы должны поочередно применять их к изображению и суммировать ответы. Обычно используется только 3-6 ориентаций, но также часто приходится фильтровать в нескольких масштабах (разные размеры ядра), поэтому в этом контексте используется большее количество фильтров.

Вы можете сделать все с помощью кода следующим образом:

img = im2double(rgb2gray(img)); %# 
[m,n] = size(img); %# Store the original size.

%# It is best if the filter size is odd, so it has a discrete center.
R = 127; C = 127;

%# The minimum amount of padding is just "one side" of the filter.
%# We add 1 if the image size is odd.
%# This assumes the filter size is odd.
pR = (R-1)/2;
pC = (C-1)/2;
if rem(m,2) ~= 0; pR = pR + 1; end;
if rem(n,2) ~= 0; pC = pC + 1; end;
img = padarray(img,[pR pC],'pre'); %# Pad image to handle circular convolution.

GW = {}; %# First, construct the filter bank.
for v = 0 : 4
    for u = 1 : 8
        GW =  [GW {GaborWavelet(R, C, Kmax, f, u, v, Delt2)}];
    end
end

%# Pad all the filters to size of padded image.
%# We made sure padsize will only be even, so we can divide by 2.
padsize = size(img) - [R C];
GW = cellfun( ...
        @(x) padarray(x,padsize/2), ...
        GW, ...
        'UniformOutput',false);

imgFFT = fft2(img); %# Pre-calculate image FFT.

for i=1:length(GW)
    filter = fft2( ifftshift( GW{i} ) ); %# See Numerical Recipes.
    imgfilt{i} = ifft2( imgFFT .* filter ); %# Apply Convolution Theorem.
end

%# Sum the responses to each filter. Do it in the above loop to save some space.
imgS = zeros(m,n);
for i=1:length(imgfilt)
    imgS = imgS + imgfilt{i}(pR+1:end,pC+1:end); %# Just use the valid part.
end

%# Look at the result.
imagesc(abs(imgS));

Имейте в виду, что это, по сути, минимальная реализация. Вы можете выбрать прокладку с репликами границы вместо нулей, применить функцию окна к изображению или увеличить размер пэда, чтобы получить частотное разрешение. Каждое из них является стандартным дополнением к описанной выше технике и должно быть тривиальным для исследования через Google и Wikipedia. Также обратите внимание, что я не добавил никаких базовых оптимизаций MATLAB, таких как предварительное выделение и т.д.

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

Ответ 2

Чтобы "применить" вейвлет к изображению, вы, как правило, берете скалярный продукт вейвлета и изображение, чтобы получить одно число, величина которого отражает то, насколько релевантен этот вейвлет для изображения. Если у вас есть полный набор всплесков (называемый "ортонормированным основанием" ) для изображения из 128 строк и 128 столбцов, у вас будет 128 * 128 = 16 384 разных всплеска. У вас здесь всего 40, но вы работаете с тем, что у вас есть.

Чтобы получить коэффициент вейвлетов, вы можете взять изображение, скажем это:

t = linspace(-6*pi,6*pi,128);
myImg = sin(t)'*cos(t) + sin(t/3)'*cos(t/3);

и возьмем скалярное произведение этого и одного из базисных векторов GW следующим образом:

myCoef = GW(:)'*myImg(:);

Мне нравится складывать все мои вейвлеты в матрицу GW_ALL, где каждая строка является одним из 32-мегапиксельных (:) 'вейвлетов, которые у вас есть, а затем вычисляет все вейвлет-коэффициенты сразу, записывая

waveletCoefficients = GW_ALL*myImg(:);

Если вы построите их со стеблем (abs (waveletCoefficients)), вы заметите, что некоторые из них больше других. Большие значения - это те, которые соответствуют изображению.

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

newImage = real(GW_ALL'*waveletCoefficients);

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

Я добавил к вашему коду (ниже), чтобы получить следующие результаты:

Если модификации:

% function gaborTest()

close all;
clear all;
clc;

% Parameter Setting
R = 128;
C = 128;
Kmax = pi / 2;
f = sqrt( 2 );
Delt = 2 * pi;
Delt2 = Delt * Delt;

% GW_ALL = nan(32, C*R);

% Show the Gabor Wavelets
for v = 0 : 4
    for u = 1 : 8
        GW = GaborWavelet ( R, C, Kmax, f, u, v, Delt2 ); % Create the Gabor wavelets
          figure( 2 );
         subplot( 5, 8, v * 8 + u ),imshow ( real( GW ) ,[]); % Show the real part of Gabor wavelets

         GW_ALL( v*8+u, :) = GW(:);

    end

    figure ( 3 );
     subplot( 1, 5, v + 1 ),imshow ( abs( GW ),[]); % Show the magnitude of Gabor wavelets

end

%% Create an Image:
t = linspace(-6*pi,6*pi,128);
myImg = sin(t)'*cos(t) + sin(t/3)'*cos(t/3);
figure(3333);
clf
subplot(1,3,1);
imagesc(myImg);
title('My Image');
axis image

%% Get the coefficients of the wavelets and plot:
waveletCoefficients = GW_ALL*myImg(:);

subplot(1,3,2);
stem(abs(waveletCoefficients));
title('Wavelet Coefficients')

%% Try and recreate the image from just a few wavelets.
% (we would need C*R wavelets to recreate perfectly)

subplot(1,3,3);
imagesc(reshape(real(GW_ALL'*waveletCoefficients),128,128))
title('My Image Reproduced from Wavelets');
axis image

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

 figure(34);
 subplot(5,8, v * 8 + u );
 imagesc(abs(ifft2((fft2(GW).*fft2(myImg)))));
 axis off