Double sum of a series

I am trying to implement the 2D convolution formula in Matlab without using conv2() built-in function. The formula is this:

where z is an [N1xN2] matrix, x is an [M1xM2] matrix and y is a [P1xP2] matrix.

N1=M1+P1-1 and N2=M2+P2-1

I did it with four for-loops but it is really slow for large matrices. Is there an faster way to do it? Here is my code:

[M1,M2] = size(x);
[P1,P2] = size(y);
 N1 = M1+P1-1;
 N2 = M2+P2-1;
 z = zeros([N1 N2]);
 y = padarray(y,[N1-P1 N2-P2],'post');

   for n1=1:N1
        for n2=1:N2
             sum1 = 0;
             for k1=1:M1
                 for k2=1:M2
                     if(n1-k1+1>0 && n2-k2+1>0)
                         sum1 = sum1 + x(k1,k2)*y(n1-k1+1,n2-k2+1);
                     end
                 end
             end
             z(n1,n2) = sum1;
         end
     end

Answer 1

[Andrei Bobrov]

function out = conv2_without_conv2(xinp,yinp)
    sx = size(xinp);
    sy = size(yinp);
    x = xinp(end:-1:1,end:-1:1);
    k = sx-1;
    so = sy + k*2;
    y = zeros(so);
    y(sx(1):end-sx(1)+1,sx(2):end-sx(2)+1) = yinp;
    i0 = reshape(1:numel(y),size(y));
    ii = i0(1:end-k(1),1:end-k(2));
    i1 = bsxfun(@plus,(0:sx(1)-1)',(0:sx(2)-1)*so(1));
    out = reshape(y(bsxfun(@plus,ii(:),i1(:)'))*x(:),sx+sy-1);
end

[xrysa]

This works perfectly! Thank you very much. I wanted to ask you about the computational complexity of your code. The one I wrote with the 4 for-loops has a complexity of O(N^4). I suppose your code should have a smaller complexity as it runs significantly faster but how do I compute it?

[xrysa]

I encounter one problem when the matrices are quite big (e.g 512 by 512). It says there is an "Error using bsxfun.Out of memory.". Is there a way to overcome this problem?

Answer 2

[Alessandro Masullo]

Why do you want to implement you own convolution when Matlab already has a very fast function for that?

For loops are slow.

Nested for loops are very slow.

Nested for loops in nested for loops are deadly slow.

If you really want to implement your own function for the convolution, and if you really want it to be fast, you need to code it in a mex file. The question is, do you really need your own implementation of conv2?

[xrysa]

Yes, I need to do it like this and then compare it to conv2 function Matlab has. I am not familiar however with coding in mex. Isn't there a way to vectorize some of the loops or the if statement to make it faster? Thanks!