beta = 1E-4; % filter parameter

% read image, prepare grid description
img = imread ( 'filter2d.pgm' );
[ N, M ] = size ( img );
h = 1 / ( N - 1 );
h2 = h * h;

% reset
Li = zeros ( 5 * N * N, 1 ); Lj = zeros ( 5 * N * N, 1 ); Lv = zeros ( 5 * N * N, 1 );
r = zeros ( N * N, 1 ); 
l = 0;

% fill matrix and right hand side
for i = 1 : N
  for j = 1 : N
    k = i + ( j - 1 ) * N;
    if i == 1 | i == N | j == 1 | j == N
      % point on boundary
      l = l + 1; Li ( l ) = k; Lj ( l ) = k; Lv ( l ) = 1; 
      r ( k ) = img ( i, j );
    else
      % point in interior
      l = l + 1; Li ( l ) = k; Lj ( l ) = k; Lv ( l ) = 4 * beta / h2 + 1; 
      l = l + 1; Li ( l ) = k; Lj ( l ) = k - 1; Lv ( l ) = - beta / h2;
      l = l + 1; Li ( l ) = k; Lj ( l ) = k + 1; Lv ( l ) = - beta / h2;
      l = l + 1; Li ( l ) = k; Lj ( l ) = k - N; Lv ( l ) = - beta / h2;
      l = l + 1; Li ( l ) = k; Lj ( l ) = k + N; Lv ( l ) = - beta / h2;
      r ( k ) = img ( i, j );
    end
  end
end
Li ( l + 1 : 5 * N * N ) = [];
Lj ( l + 1 : 5 * N * N ) = [];
Lv ( l + 1 : 5 * N * N ) = [];
L = sparse ( Li, Lj, Lv );

% solve and draw
u = L \ r;
res = img;
res ( : ) = u;
subplot ( 121 );
imagesc ( img );
set( gcf, 'Unit', 'pixels', 'Position', [ 10 100 2 * N + 30 N + 20 ] );
set( gca, 'Unit', 'pixels', 'Position', [ 10 10 N N ] );
axis off;
subplot ( 122 );
imagesc ( res ); 
set( gca, 'Unit', 'pixels', 'Position', [ 20 + N 10 N N ] );
axis off;
colormap ( gray );