在 matlab 中处理可分离的 Gabor 滤波器

Question

如果过滤器 g 可以表示为两个向量 grow 和 gcol 的乘积，则它称为可分离过滤器。使用一维过滤器将二维过滤器的计算复杂度从O(M^2 N^2) 降低到O(2M N^2)，其中 M 和 N 分别是过滤器掩码和图像的宽度（和高度）。

在this stackoverflow link 中，我编写了空间域中的 Gabor 滤波器方程，然后编写了一个 matlab 代码，用于创建 64 个 gabor 特征。

根据可分离过滤器的定义，Gabor 过滤器平行于图像轴 - theta = k*pi/2 where k=0,1,2,etc.。所以如果 theta=pi/2 ==> this stackoverflow link 中的等式可以改写为：

上面的等式是从this article中提取的。

注意： theta 可以扩展到等于k*pi/4. 通过比较this stackoverflow link 中的等式，我们可以认为f= 1 / lambda。

通过更改我之前在this stackoverflow link 中的代码，我编写了一个 matlab 代码，通过使用上面的等式使 Gabor 滤波器可分离，但我确信下面的代码不正确，尤其是当我初始化 @ 987654334@ 和 glp 方程。这就是为什么我需要你的帮助。非常感谢您的帮助。

现在让我们展示我的代码：

    function [fSiz,filters1,filters2,c1OL,numSimpleFilters] = init_gabor(rot, RF_siz)

    image=imread('xxx.jpg');
    image_gray=rgb2gray(image);
    image_gray=imresize(image_gray, [100 100]);
    image_double=double(image_gray);

     rot = [0 45 90 135]; % we have four orientations
                    RF_siz    = [7:2:37]; %we get 16 scales (7x7 to 37x37 in steps of two pixels)
                    minFS     = 7; % the minimum receptive field
                    maxFS     = 37; % the maximum receptive field
                    sigma  = 0.0036*RF_siz.^2 + 0.35*RF_siz + 0.18; %define the equation of effective width
                    lambda = sigma/0.8; % it the equation of wavelength (lambda)
                    G      = 0.3;   % spatial aspect ratio: 0.23 < gamma < 0.92


                    numFilterSizes   = length(RF_siz); % we get 16

                    numSimpleFilters = length(rot); % we get 4

                    numFilters       = numFilterSizes*numSimpleFilters; % we get 16x4 = 64 filters

                    fSiz             = zeros(numFilters,1); % It is a vector of size numFilters where each cell contains the size of the filter (7,7,7,7,9,9,9,9,11,11,11,11,......,37,37,37,37)
filters1          = zeros(max(RF_siz),numFilters);
filters2          = zeros(numFilters,max(RF_siz));

for k = 1:numFilterSizes  
    for r = 1:numSimpleFilters
        theta     = rot(r)*pi/180;
        filtSize  = RF_siz(k);
        center    = ceil(filtSize/2);
        filtSizeL = center-1;
        filtSizeR = filtSize-filtSizeL-1;
        sigmaq    = sigma(k)^2;

        for x = -filtSizeL:filtSizeR

                    fx = exp(-(x^2)/(2*sigmaq))*cos(2*pi*x/lambda(k));
                       f1(x+center,1) = fx;
        end
                      for y = -filtSizeL:filtSizeR
                    gy = exp(-(y^2)/(2*sigmaq));
                    f2(1,y+center) = gy;
                      end



        f1 = f1 - mean(mean(f1));
        f1 = f1 ./ sqrt(sum(sum(f1.^2)));
         f2 = f2 - mean(mean(f2));
        f2 = f2 ./ sqrt(sum(sum(f2.^2)));
        p = numSimpleFilters*(k-1) + r;
        filters1(1:filtSize,p)=f1;
        filters2(p,1:filtSize)=f2;

        convv1=imfilter(image_double,  filters1(1:filtSize,p),'conv');
        convv2=imfilter(double(convv1),  filters2(p,1:filtSize),'conv');
        figure
        imagesc(convv2);
        colormap(gray);

    end
end

Answer 1

最终代码为：

function [fSiz,filters1,filters2,c1OL,numSimpleFilters] = init_gabor(rot, RF_siz)

            image=imread('xxx.jpg');
            image_gray=rgb2gray(image);
            image_gray=imresize(image_gray, [100 100]);
            image_double=double(image_gray);

             rot = [0 45 90 135];
                            RF_siz    = [7:2:37]; 
                            minFS     = 7; 
                            maxFS     = 37; 
                            sigma  = 0.0036*RF_siz.^2 + 0.35*RF_siz + 0.18; 
                            lambda = sigma/0.8; 
                            G      = 0.3;   


                            numFilterSizes   = length(RF_siz); 

                            numSimpleFilters = length(rot); 

                            numFilters       = numFilterSizes*numSimpleFilters; 

                            fSiz             = zeros(numFilters,1); 
                            filters1          = zeros(max(RF_siz),numFilters);
                            filters2          = zeros(numFilters,max(RF_siz));

        for k = 1:numFilterSizes  
            for r = 1:numSimpleFilters
                theta     = rot(r)*pi/180;
                filtSize  = RF_siz(k);
                center    = ceil(filtSize/2);
                filtSizeL = center-1;
                filtSizeR = filtSize-filtSizeL-1;
                sigmaq    = sigma(k)^2;

                for x = -filtSizeL:filtSizeR

                            fx = exp(-(x^2)/(2*sigmaq))*exp(sqrt(-1)*x*cos(theta));
                               f1(1, x+center) = fx;
                end
                              for y = -filtSizeL:filtSizeR
                            gy=exp(-(y^2)/(2*sigmaq))*exp(sqrt(-1)*y*sin(theta));
                            f2(y+center,1) = gy;
                              end


         f1 = f1 - mean(mean(f1));
            f1 = f1 ./ sqrt(sum(sum(f1.^2)));
             f2 = f2 - mean(mean(f2));
            f2 = f2 ./ sqrt(sum(sum(f2.^2)));
            p = numSimpleFilters*(k-1) + r;

            filters1(1:filtSize,p)=f1;
            filters2(p,1:filtSize)=f2;

            convv1=imfilter(image_double,  filters1(1:filtSize,p),'conv');
            convv2=imfilter(double(convv1),  filters2(p,1:filtSize),'conv');

            figure
            imagesc(imag(convv2));
            colormap(gray);

        end
    end

Answer 2

如果您以前版本的 Gabor 过滤器代码也是正确的，我认为代码是正确的。唯一的问题是，如果theta = k * pi/4;，您的公式here 应分隔为：

fx = exp(-(x^2)/(2*sigmaq))*cos(2*pi*x/lambda(k));
gy = exp(-(G^2 * y^2)/(2*sigmaq));

为了保持一致，你可以使用

f1(1,x+center) = fx;
f2(y+center,1) = gy;

或保持 f1 和 f2 不变，但之后转置您的 filters1 和 filters2。 在我看来，其他一切都很好。

编辑

根据您的论文，我上面的答案适用于theta = k * pi/4;，还有其他角度，

x = i*cos(theta) - j*sin(theta);
y = i*sin(theta) + j*cos(theta);
fx = exp(-(x^2)/(2*sigmaq))*exp(sqrt(-1)*x*cos(theta));
gy = exp(-(G^2 * y^2)/(2*sigmaq))*exp(sqrt(-1)*y*sin(theta));