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 

正如你所看到的，外兩個環的滑動窗口，而剩下的那些內部是求和的變量。

現在，我向你請求我的親愛的計算器用戶是：你能幫我改進這些非常討厭的循環，其花費的時間超過它的份額，並把它變成向量化的循環？ 並且這種改善會有重大改變嗎？

謝謝。

Answer 1

這四個嵌套循環基本上是以滑動鄰域風格處理圖像中的每個像素。我立即想到了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,[]) 

，那麼你可能需要對矩陣進行縮放/歸一化，並應用閾值...

Answer 2

這可能不是你要求的，但似乎（乍一看）總和的順序是獨立的，而不是{x，y，l，k}你可以去{l，k，的x，y}。這樣做將允許您通過將內核保留在最外層循環中來更少地評估內核。