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

Posted 2023-03-15

【中文标题】如何将 Gabor 小波应用于图像？

【英文标题】：How to apply Gabor wavelets to an image?

【发布时间】：2012-02-18 15:54:18

【问题描述】：

如何在图像上应用这些 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

编辑：这是我图片的尺寸

【问题讨论】：

如果尺寸给您带来麻烦，您可以考虑设置 R=size(img,1); C=size(img,2) 并使用 img=double(rgb2gray(img)) 转换为灰度并对其进行计算。

您使用的是彩色图像，它是一个 3D 数组。您应该阅读有关基本image representation 的教程。

我确实将其转换为灰度图像

您好，如何应用逆Gabor小波公式来重建原始图像？任何帮助将不胜感激

克里斯蒂娜，你的问题是另一个问题，应该单独发布。不过，请随时在评论中链接到您的新问题，因为它在这里非常相关。

【参考方案1】：

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

您可以使用Convolution Theorem 对图像与图像进行卷积，方法是对图像和过滤器的傅里叶变换的逐元素乘积进行傅里叶逆变换。这是基本公式：

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

通常，FFT 卷积对于大小 > 20 的内核更为出色。有关详细信息，我推荐 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( GWi ) ); %# See Numerical Recipes.
    imgfilti = 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 + imgfilti(pR+1:end,pC+1:end); %# Just use the valid part.
end

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

请记住，这基本上是最低限度的实现。您可以选择使用复制的边界而不是零来填充，对图像应用窗口函数或使填充尺寸更大以获得频率分辨率。这些中的每一个都是对我上面概述的技术的标准增强，通过Google 和Wikipedia 进行研究应该是微不足道的。另请注意，我没有添加任何基本的 MATLAB 优化，例如预分配等。

最后一点，如果您的过滤器总是比图像小很多，您可能希望跳过图像填充（即使用第一个代码示例）。这是因为，向图像添加零会在填充开始处创建人工边缘特征。如果过滤器很小，则循环卷积的环绕不会导致问题，因为只会涉及 过滤器 填充中的零。但是一旦过滤器足够大，环绕效应就会变得严重。如果您必须使用大型滤镜，您可能需要使用更复杂的填充方案，或者裁剪图像的边缘。

【讨论】：

我把提到的代码放在哪里：[m,n] = size(img); [mf,nf] = 尺寸（GW）； GW = padarray(GW,[n-nf m-mf]/2); GW = ifftshift(GW); imgf = ifft2( fft2(img) .* fft2(img) );

该代码只是给出基本公式。我进行了编辑以使其更加清晰，并添加了更多信息。第二个代码块应该覆盖它。

我在 padarray 中遇到错误？？？使用 ==> iptcheckinput 函数时出错 PADARRAY 期望其第二个输入 PADSIZE 为整数值。 ==> padarray>ParseInputs 中的错误 233 iptcheckinput(padSize, 'double', 'real' 'vector' 'nonnan' 'nonnegative' ... ==> padarray 中的错误在 65 [a, method, padSize , padVal, direction] = ParseInputs(varargin:); ==> @(x)padarray(x,[m,n]/2) 错误 ==> GaborExample at 28 GW = cellfun(@(x )padarray(x,[mn]/2),GW,'UniformOutput',false);

???索引超出矩阵维度。 ==> GaborExample 中的错误在 32 GW = [GWGaborWavelet(R,C,Kmax,f,u,v,Delt2)];我不知道为什么我会收到这个错误

我已将颜色转换代码添加到我的第一个代码示例中。我强烈建议您阅读有关 image representation 和操作的 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 个 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 小波（如前所述）不是正交的（reference），并且更有可能用于使用 reve_etrange 描述的卷积进行特征提取。在这种情况下，您可能会考虑将其添加到您的内部循环中：

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

【讨论】：

@vini 你想做什么？我做了一些修改，但我倾向于同意 reve 的观点，即您可能正在寻找卷积。查看参考disp.ee.ntu.edu.tw/~pujols/…

是的，这就是我真正想要的！我正在实施一篇论文

Gabor 滤波器和 Gabor 小波有什么区别？ Gabor 小波是否更优化？请感谢您的回答:) 提前致谢！

@Christina 我也想知道。我猜他们指的是同一个东西，即 40 名成员的过滤器组。还有一个问题。为什么我们取内积而不是卷积？还是它们是一样的？

小波被用作过滤器，因此在上下文中 Gabor wavelet == Gabor 过滤器。内积与卷积不同。使用内积会给出一个数字，即图像与过滤器的相似度，可用于分类。卷积给出了一个新的图像，表示与过滤器在每个点的相似性。这对于边缘检测等很有用。