Geek asks and answers

MATLAB double summation and vectorized loops

Here is my attempt to implement this beautiful formula.

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

As you can see, the outer two loops are for the sliding window, and the remaining inner ones are for summation variables.

Now, my request is for you, my favorite stackoverflow users: Can you help me improve these very unpleasant for loops that take more than its share of time and turn them into vectorized loops? And will this improvement be a significant change?

Thank.

+3
source share
2 answers

. 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

( , , ):

%# 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,[])

/ ...

+1

, , ( ), {x, y, l, k} {l, k, x, }. , .

+4

More articles:


All Articles