Skip to content

Instantly share code, notes, and snippets.

@d1manson

d1manson/so_answer.m

Created May 2, 2016 14:47
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save d1manson/37ac0968c4a79ab62ad41cc29a36dbbf to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save d1manson/37ac0968c4a79ab62ad41cc29a36dbbf to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. Daniel Manson created this gist May 2, 2016.
    27 changes: 27 additions & 0 deletions so_answer.m
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,27 @@
    % sample data
    s1 = 8192; s2 = 200;
    img_a = rand(s1, s2);
    img_b = rand(s1, s2);
    r = 2;

    % and the calculation itself
    img_diff = img_a - img_b;
    kernel = bsxfun(@plus, (-s1:s1).^2', (-s2:s2).^2);
    kernel = 1/(2/pi/r^2) * exp(-kernel/ (2*r*2));
    g = conv2(img_diff, kernel, 'same'); % see below for faster, fft based convolution
    res = g(:)' * img_diff(:);


    %%
    function c = conv2fft(X, Y)
    % ignoring small floating-point differences, this is equivalent
    % to the inbuilt Matlab conv2(X, Y, 'same').
    % https://en.wikipedia.org/wiki/Convolution_theorem
    % http://stackoverflow.com/a/21757317/2399799
    X1 = [X zeros(size(X,1), size(Y,2)-1);
    zeros(size(Y,1)-1, size(X,2)+size(Y,2)-1)];
    Y1 = zeros(size(X1));
    Y1(1:size(Y,1), 1:size(Y,2)) = Y;
    c = ifft2(fft2(X1).*fft2(Y1));
    c = c(size(X,1)+1:size(X,1)+size(X,1), size(X,2)+1:size(X,2)+size(X,2));
    end