function Kf = my_conv2(f,k)
% function Kf = my_conv2(f,k)
% convolution of image f with a kernel k
% reflection boundary conditions are used
% returns same size image as f
% note, order of input is important
% no error checking is built in

[m,n] = size(k);
[M,N] = size(f);

rx = (m-1)/2; %radius of kernel
ry = (n-1)/2;
p = M+2*rx; %size of padded matrix
q = N+2*ry; %size of padded matrix
fp = zeros(p,q);

%reflect boundary
fp(rx+1:M+rx,ry+1:N+ry) = f;
fp(rx+1:M+rx,1:ry) = f(1:M,ry:-1:1); 
fp(rx+1:M+rx,q-ry+1:q) = f(1:M,N:-1:N-ry+1);
fp(1:rx,:) = fp(2*rx:-1:rx+1,:); 
fp(p-rx+1:p,:) = fp(p-rx:-1:p-2*rx+1,:);

Kf = conv2(fp,k,'same');

Kf = Kf(rx+1:p-rx,ry+1:q-ry); % extract non-padded region

end

