function [result]=convolution(image,kernel)
%USAGE: convolution(image,kernel)
% Computes convolution of two matrices. Image can be 2D or 3D array and kernel must be 2D array.
%EXAMPLE:
% Gaussian blur convolution(imread('im.png'),(1/(16)) * [1 2 1;2 4 2;1 2 1]) 
% Laplace edge detector convolution(imread('im.png'),[1 1 1;1 -8 1;1 1 1])

imtype=class(image);

switch imtype
  case {'uint8','uint16','uint32','uint64'}  
  otherwise
    error('Image must be integer type.');
end

kern=struct('data',rot90(rot90(kernel(:,:,1))),'size',size(kernel));
if (length(kern.size)>2)
  error('Convolution kernel must be 2D array.');
end

im=struct('data',image,'size',size(image));
g=round(kern.size(1:2)/2)+1;

if (length(im.size)>2) 
  i=im.size(3);
else 
  i=1;
end

result=zeros(im.size,imtype);

for j=1:i
  ret=zeros(im.size(1:2)+(kern.size(1:2)-1),class(kernel));
  ret(g(1):g(1)+im.size(1)-1,g(2):g(2)+im.size(2)-1) = im.data(:,:,j);

  for m = [1:im.size(1)]
    for n = [1:im.size(2)]
      result(m,n,j) = cast(sum(sum(ret(m:m+kern.size(1)-1,n:n+kern.size(2)-1).*kern.data)),imtype);
    end
  end
end
