MATLAB 和矢量化循环中的双重求和

Question

这是我实现这个可爱公式的尝试。

http://dl.dropbox.com/u/7348856/Picture1.png

%WIGNER Computes Wigner-Distribution on an image (difference of two images).
function[wd] = wigner(difference)
%Image size
[M, N, ~] = size(difference);
%Window size (5 x 5)
Md = 5;
Nd = 5;
%Fourier Transform
F = fft2(difference);
%Initializing the wigner picture
wd = zeros(M, N, 'uint8');
lambda =0.02;
value = (4/(Md*Nd));
for x = 1+floor(Md/2):M - floor(Md/2)
    for y = 1+floor(Nd/2):N - floor(Nd/2)
        for l = -floor(Nd/2) : floor(Nd/2)
            for k = -floor(Md/2) : floor(Md/2)
                kernel = exp(-lambda * norm(k,l));
                kernel = kernel * value;
                theta = 4 * pi * ((real(F(x, y)) * (k/M) )+ (imag(F(x, y)) * (l/N)));
                wd(x, y) = (wd(x, y)) + (cos(theta) * difference(x + k, y + l) * difference(x - k, y - l) * (kernel));
            end
        end
    end
end
end

如您所见，外面的两个循环用于滑动窗口，而其余的内部循环用于求和的变量。

现在，亲爱的 stackoverflow 用户，我对您的请求是：您能帮我改进这些非常讨厌的 for 循环，这些循环花费的时间超过其份额，并将其变成矢量化循环吗？ 这种改进会带来重大变化吗？

谢谢。

Answer 1

这可能不是您要问的，但似乎（乍一看）求和的顺序是独立的，您可以用 {l,k,x 代替 {x,y,l,k} ,y}。这样做可以让您通过将内核保持在最外层循环中来减少对内核的评估次数。

Answer 2

这四个嵌套循环基本上是以滑动邻域样式处理图像中的每个像素。我立刻想到了NLFILTER 和IM2COL 函数。

这是我对代码进行矢量化的尝试。请注意，我尚未对其进行彻底测试，或将性能与基于循环的解决方案进行比较：

function WD = wigner(D, Md, Nd, lambda)
    %# window size and lambda
    if nargin<2, Md = 5; end
    if nargin<3, Nd = 5; end
    if nargin<4, lambda = 5; end

    %# image size
    [M,N,~] = size(D);

    %# kernel = exp(-lambda*norm([k,l])
    [K,L] = meshgrid(-floor(Md/2):floor(Md/2), -floor(Nd/2):floor(Nd/2));
    K = K(:); L = L(:);
    kernel = exp(-lambda .* sqrt(K.^2+L.^2));

    %# frequency-domain part
    F = fft2(D);

    %# f(x+k,y+l) * f(x-k,y-l) * kernel
    C = im2col(D, [Md Nd], 'sliding');
    X1 = bsxfun(@times, C .* flipud(C), kernel);

    %# cos(theta)
    C = im2col(F, [Md Nd], 'sliding');
    C = C(round(Md*Nd/2),:);    %# take center pixels
    theta = bsxfun(@times, real(C), K/M) + bsxfun(@times, imag(C), L/N);
    X2 = cos(4*pi*theta);

    %# combine both parts for each sliding-neighborhood
    WD = col2im(sum(X1.*X2,1), [Md Nd], size(F), 'sliding') .* (4/(M*N));

    %# pad array with zeros to be of same size as input image
    WD = padarray(WD, ([Md Nd]-1)./2, 0, 'both');
end

对于它的价值，这里是基于循环的版本，具有@Laurbert515 建议的改进：

function WD = wigner_loop(D, Md, Nd, lambda)
    %# window size and lambda
    if nargin<2, Md = 5; end
    if nargin<3, Nd = 5; end
    if nargin<4, lambda = 5; end

    %# image size
    [M,N,~] = size(D);

    %# frequency-domain part
    F = fft2(D);

    WD = zeros([M,N]);
    for l = -floor(Nd/2):floor(Nd/2)
        for k = -floor(Md/2):floor(Md/2)
            %# kernel = exp(-lambda*norm([k,l])
            kernel = exp(-lambda * norm([k,l]));

            for x = (1+floor(Md/2)):(M-floor(Md/2))
                for y = (1+floor(Nd/2)):(N-floor(Nd/2))
                    %# cos(theta)
                    theta = 4 * pi * ( real(F(x,y))*k/M + imag(F(x,y))*l/N );

                    %# f(x+k,y+l) * f(x-k,y-l)* kernel
                    WD(x,y) = WD(x,y) + ( cos(theta) * D(x+k,y+l) * D(x-k,y-l) * kernel );
                end
            end
        end
    end
    WD = WD * ( 4/(M*N) );
end

以及我如何测试它（基于我从paper 你previously 链接到的内容）：

%# difference between two consecutive frames
A = imread('AT3_1m4_02.tif');
B = imread('AT3_1m4_03.tif');
D = imsubtract(A,B);
%#D = rgb2gray(D);
D = im2double(D);

%# apply Wigner-Distribution
tic, WD1 = wigner(D); toc
tic, WD2 = wigner_loop(D); toc
figure(1), imshow(WD1,[])
figure(2), imshow(WD2,[])

然后您可能需要对矩阵进行缩放/归一化，并应用阈值...