3

これが、この素敵な式を実装するための私の試みです。

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

ご覧のとおり、外側の2つのループはスライディングウィンドウ用で、残りの内側のループは合計の変数用です。

さて、私の最愛のstackoverflowユーザーへの私のリクエストは、次のとおりです。時間のシェア以上の時間がかかるこれらの非常に厄介なforループを改善し、ベクトル化されたループに変換するのを手伝ってくれませんか？そして、その改善は大きな変化になるのでしょうか？

ありがとうございました。

4

2 に答える 2

4

これはあなたが求めているものではないかもしれませんが、（一見すると）合計の順序は独立しており、{x、y、l、k}の代わりに{l、k、x、y}に行くことができるようです。 。これを行うと、カーネルを最も外側のループに保持することで、カーネルを評価する回数を減らすことができます。

于 2012-05-25T17:42:30.147 に答える
1

これらの4つのネストされたループは、基本的に画像内の各ピクセルをスライディングネイバーフッドスタイルで処理しています。私はすぐに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,[])

次に、マトリックスをスケーリング/正規化し、しきい値を適用する必要がある場合があります...

于 2012-05-27T00:59:07.710 に答える