如何将Gabor小波应用于图像？

Question

我如何在图像上应用这些 Gabor 滤波小波？

close all;
clear all;
clc;

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

% 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

    end

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

end



function GW = GaborWavelet (R, C, Kmax, f, u, v, Delt2)


k = ( Kmax / ( f ^ v ) ) * exp( 1i * u * pi / 8 );% Wave Vector

kn2 = ( abs( k ) ) ^ 2;

GW = zeros ( R , C );

for m = -R/2 + 1 : R/2

    for n = -C/2 + 1 : C/2

        GW(m+R/2,n+C/2) = ( kn2 / Delt2 ) * exp( -0.5 * kn2 * ( m ^ 2 + n ^ 2 ) / Delt2) * ( exp( 1i * ( real( k ) * m + imag ( k ) * n ) ) - exp ( -0.5 * Delt2 ) );

    end

end

编辑：这是我的图片尺寸

- vini

---

如果尺寸让您感到困扰，您可以考虑设置R=size(img,1); C=size(img,2)，并使用img=double(rgb2gray(img))将其转换为灰度图像，并在此基础上进行计算。 - Steve

您正在使用彩色图像，它是一个三维数组。您应该阅读有关基本图像表示的教程。 - reve_etrange

我已将其转换为灰度图像。 - vini

你好，如何应用逆 Gabor 小波公式来重建原始图像？非常感谢任何帮助。 - Christina

Christina，你的问题与这个不同，应该单独发布。但是，如果你觉得它在这里非常相关，可以在评论中链接到你的新问题。 - reve_etrange

@Christina 提出问题会帮助你解决问题。 - vini

Answer 1

Gabor滤波器的典型用途是在多个方向上计算滤波器响应，例如用于边缘检测。

您可以使用卷积定理通过对图像和滤波器的傅里叶变换的逐元素乘积的傅里叶反变换来将滤波器与图像进行卷积。以下是基本公式：

%# 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) .* GW );

通常情况下，对于大于20的内核大小，FFT卷积优于其他方法。如需详细了解，请参考《Numerical Recipes in C》，其中有关于该方法及其注意事项的良好、语言无关的描述。

您的内核已经很大了，但使用FFT方法后，它们可以达到与图像一样的大小，因为它们会被填充到那个尺寸。由于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优化，例如预分配等。

最后，如果您的滤波器始终比图像小得多，则可能要跳过图像填充（即使用第一个代码示例）。这是因为向图像添加零会创建人工边缘特征，这是填充开始的地方。如果滤波器很小，来自循环卷积的环绕不会引起问题，因为只有在滤波器的填充中的零会参与其中。但是，一旦滤波器足够大，环绕效应就会变得严重。如果必须使用大型滤波器，则可能必须使用更复杂的填充方案或裁剪图像的边缘。

- reve_etrange

---

我应该把这段代码放在哪里： [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) ); - vini

我在padarray中遇到了一个错误 ??? Error using ==> iptcheckinput Function PADARRAY expected its second input, PADSIZE, to be integer-valued.Error in ==> padarray>ParseInputs at 233 iptcheckinput(padSize, {'double'}, {'real' 'vector' 'nonnan' 'nonnegative' ...Error in ==> padarray at 65 [a, method, padSize, padVal, direction] = ParseInputs(varargin{:});Error in ==> @(x)padarray(x,[m,n]/2)Error in ==> GaborExample at 28 GW = cellfun(@(x)padarray(x,[m n]/2),GW,'UniformOutput',false); - vini

??? 索引超出矩阵维度。错误位于 ==> GaborExample 第32行 GW = [GW{GaborWavelet(R,C,Kmax,f,u,v,Delt2)}]; 我不知道为什么会出现这个错误。 - vini

我正在尝试，但问题一再出现。 - vini

我已经将颜色转换代码添加到我的第一个代码示例中。我强烈建议您阅读MATLAB关于图像表示和操作的文档。 - reve_etrange

Answer 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个GW(:)小波之一，然后通过编写代码一次性计算所有小波系数。

waveletCoefficients = GW_ALL*myImg(:);

如果您使用stem(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

这种方法是提取小波系数和重建图像的基础。正如参考文献所述，Gabor小波不是正交的，更适合使用卷积进行特征提取。在这种情况下，你可以考虑将其添加到你的内部循环中：

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

- Steve

---

@vini 你想做什么？我进行了一些小的编辑，但我倾向于同意reve的看法，你可能在寻找卷积。请查看参考文献http://disp.ee.ntu.edu.tw/~pujols/Gabor%20wavelet%20transform%20and%20its%20application.pdf - Steve

是的，那正是我正在寻找的！我正在实现一篇论文。 - vini

Gabor滤波器和Gabor小波之间有什么区别？ Gabor小波是否更优？请回答，谢谢！ - Christina

@Christina 我也在想同样的问题。我猜他们指的是同一个东西，那个由40个成员组成的滤波器组。还有一个问题，为什么我们要取内积而不是卷积？或者它们是一样的吗？ - jeff

小波被用作滤波器，因此在上下文中Gabor小波== Gabor滤波器。内积与卷积不同。使用内积会得到一个单一的数字，即图像与滤波器的相似度，这可能用于分类。卷积会给出一个新的图像，指示每个点与滤波器的相似度。这对于边缘检测等操作非常有用。 - reve_etrange