Code covered by the BSD License

# Matrix Convolution with Sub-Pixel Resolution

### Tristan Ursell (view profile)

Convolve two matrices at sub-pixel resolution, using bilinear interpolation.

matout=matoverlay(mat1,mat2,x,y)
```function matout=matoverlay(mat1,mat2,x,y)
%Tristan Ursell
%Sub-pixel Resolved 2D Convolution
%March 2012
%
%matout=matoverlay(mat1,mat2,x,y);
%
%This function takes an input matrix mat1 and creates an image of the
%matrix mat2 at the position (x,y) in mat1.  If (x,y) are floats, then the
%image is a sub-pixel bilinear representatoin of mat2 at position (x,y) in
%mat1. The output matrix will have the same size at mat1, with no edge
%effects.
%
%Essentially this is performing a sparse, fully valid convolution of mat2
%and mat1 at the point (x,y) with the output size of mat1.  The point (x,y)
%uses the imaging convention for the coordinate axes.
%
%The values of (x,y) can be floats, as long as they lie within the bounds
%of mat1.  Combining this function with a for-loop and weights creates a
%fully valid 2D subpixel resolution convolution -- see Example -- in
%contrast to conv2 which is limited to pixel resolution.
%
%
%Example:
%
%N=50;
%x=1+99*rand(1,N);
%y=1+99*rand(1,N);
%
%mat1=zeros(100,100);
%
%mat2=mat2gray(fspecial('gaussian',[11,11],3));
%
%I0=zeros(size(mat1));
%ints=rand(1,N);
%for i=1:N
%   I0=I0+ints(i)*matoverlay(mat1,mat2,x(i),y(i));
%end
%
%figure;
%colormap(hot)
%subplot(1,2,1)
%imagesc(mat2)
%axis equal
%axis tight
%title('mat2')
%
%subplot(1,2,2)
%hold on
%imagesc(I0)
%plot(x,y,'bo')
%axis equal
%axis tight
%title('Sparse Convolution of mat1 and mat2')
%

%get sizes of input matrices
[sy1,sx1]=size(mat1);
[sy2,sx2]=size(mat2);

%check the values of x and y
if or(floor(x)<1,floor(y)<1)
error('The coordinate point must round to a value greater than (1,1).')
end
if or(ceil(x)>sx1,ceil(y)>sy1)
error('The coordinate point must round to a valid value.')
end

%pad mat2 with zeros
mat2b=zeros(sx2+2,sy2+2);
mat2b(2:end-1,2:end-1)=mat2;
clear mat2
mat2=mat2b;
[sy2,sx2]=size(mat2);

%find indices of mat2 center
xmid=round(sx2/2);
ymid=round(sy2/2);

%revert to pixel-resolution
%x=round(x);
%y=round(y);

%fix numerical precision offset
if floor(x)==x
x=x+1e-12;
end
if floor(y)==y
y=y+1e-12;
end

%calculate weights (starting from top left, clockwise)
w1=(x-floor(x))*(y-floor(y));
w2=(ceil(x)-x)*(y-floor(y));
w3=(ceil(x)-x)*(ceil(y)-y);
w4=(x-floor(x))*(ceil(y)-y);

%construct weighted kernel
kernel=zeros(size(mat2));
kernel(2:end,2:end)=kernel(2:end,2:end)+w1*mat2(1:end-1,1:end-1);
kernel(2:end,:)=kernel(2:end,:)+w2*mat2(1:end-1,:);
kernel=kernel+w3*mat2;
kernel(:,2:end)=kernel(:,2:end)+w4*mat2(:,1:end-1);

%create temporary new mat1 with larger mat2 sized borders.
mattemp=zeros(sy1+2*ymid+2,sx1+2*xmid+2);

%place down green's function
mattemp(floor(y)+1:floor(y)+sy2,floor(x)+1:floor(x)+sx2)=kernel;

%set output matrix
matout=mattemp(ymid+1:ymid+sy1,xmid+1:xmid+sx1);

```