clear
close all

%% Parameter
beta=1e-4;

%% Read data
display('reading file');
inputimage = imread('filter2d.pgm');
f=double(inputimage)/255;
[n,m] = size(f);
N = m*n;
h = 1/m;

%% Display input
image(f,'CDataMapping','scaled')
colormap gray
title('Input image')
axis equal
axis off
drawnow

%% Assemble matrices, column-wise node counting
display('assembling');
e=ones(N,1);
% Laplace
A=spdiags([e e -4.*e e e],[-n -1 0 1 n],N,N);
% boundary conditions
% left boundary
display('upper bc');
A(1:n,:)=0;
% right boundary
display('lower bc');
A(N-n+1:N,:)=0;
% upper boundary
display('left bc');
A(1:m:N,:)=0;
% lower boundary
display('right bc');
A(m:m:N,:)=0;

B=spdiags(e,0,N,N);
% Operator
OP=-beta/(h*h)*A + B;
rhs = f(:);

%% Show matrix
figure
spy(OP)

%% Solve
display('solving');
sol=OP\rhs;
solution=reshape(sol,n,m);

%% Plot
figure
title('Output image')
image(solution,'CDataMapping','scaled')
colormap gray
axis equal
axis off
%
imwrite(solution,'filter2d_result.pgm','pgm');