clear; clc; clf;

N = 25;
m = N+1;
h = 1/N+1; my = 1.7;
N4 = round(N/4);
ustar = [zeros(N4,1); ones(m-2*N4,1); zeros(N4,1)];
ustar = kron(ustar, ustar');
v = ustar + 0.1*randn(m,m);

M = m*m;
v = reshape(v, M, 1);

%% D erstellen:
% b0 ... Hauptdiagonale
b0 = zeros(M, 1);
for i = 1:(m-2)
    b0(i*m+1:(i*m+m), 1) = [3; 4*ones(m-2,1); 3];
end
b0(1:m, 1) = [2; 3*ones(m-2,1); 2];
b0(m^2-m+1:end, 1) = [2; 3*ones(m-2,1); 2];

% b1 ...  Nebendiagonalen +1, -1
b1 = -ones(M, 1);
for i = 0:(m-1)
    b1(i*m+1,1) = 0;
end
D = spdiags([b0, b1,[b1(2:end,1);-1],-ones(M, 1),-ones(M, 1)], ...
    [0, 1, -1, m, -m], m^2, m^2);

%% GAUSS (ohne Speichern der Multiplikatoren)
tic

A = full([(my/h^2*D + speye(m^2)), v]);
for k = 1:(M-1)                % Pivotelement
	for i=(k+1):min(M, k+m)    % Zeilenindex
		x = A(i, k)/A(k, k);   % Multiplikator
		% das x-fache der k-ten Zeile von der i-ten abziehen
        % dabei auf v nicht vergessen!
        A(i, [k+1:min(M, k+m), M+1]) =   ...
          A(i, [k+1:min(M, k+m), M+1]) - x * A(k, [k+1:min(M, k+m), M+1]);
	end
end
disp(['Matrix-Dimension: ', num2str(M), ' x ', num2str(M)]);
disp('Zeit für GAUSS:');
toc

%% Rückwärts Einsetzen
u = full([zeros(M-1, 1); A(M,M+1)/A(M,M)]);
for i = M-1:-1:1
    u(i) = (A(i, M+1) - dot(A(i, i+1:min(i+m,M)).', u(i+1:min(i+m,M)))) ...
           /A(i,i);
end

%% Probe
u2 = (my/h^2*D + speye(m^2))\v;
disp(['PROBE: norm(u-u2,1) = ', num2str(norm(u-u2,1))]);

%% Plot
u = reshape(u, m, m);
v = reshape(v, m, m);

h1 = figure(1); clf; set(h1,'Position',[10 10 900 300]);
subplot(1,3,1);
imagesc(ustar); colormap('gray'); axis image; axis off;
title('exakt');
subplot(1,3,2);
imagesc(v); colormap('gray'); axis image; axis off;
title('verrauscht');
subplot(1,3,3);
imagesc(u); colormap('gray'); axis image; axis off;
title('entrauscht');

input('Setzte mit echten Daten fort?')

clear; clc; clf;

v = imread('Pappos.bmp');
v = v(:,:,1);
v = double(v);
N = size(v,1)-1;
m = N+1;
h = 1/m; 
my = 4e-6;
N4 = round(N/4);

M = m*m;
v = reshape(v, M, 1);

%% D erstellen:
% b0 ... Hauptdiagonale
b0 = zeros(M, 1);
for i = 1:(m-2)
    b0(i*m+1:(i*m+m), 1) = [3; 4*ones(m-2,1); 3];
end
b0(1:m, 1) = [2; 3*ones(m-2,1); 2];
b0(m^2-m+1:end, 1) = [2; 3*ones(m-2,1); 2];

% b1 ...  Nebendiagonalen +1, -1
b1 = -ones(M, 1);
for i = 0:(m-1)
    b1(i*m+1,1) = 0;
end
D = spdiags([b0, b1,[b1(2:end,1);-1],-ones(m*m, 1),-ones(m*m, 1)], ...
    [0, 1, -1, m, -m], m^2, m^2);


%% Lösen (mit Backslash, weil schneller)
u = (my/h^2*D + speye(m^2)) \ v;

%% Plot
u = reshape(u, m, m);
v = reshape(v, m, m);

h1 = figure(1); clf; set(h1,'Position',[10 10 900 300]);
subplot(1,2,1);
imagesc(v); colormap('gray'); axis image; axis off;
title('original');
subplot(1,2,2);
imagesc(u); colormap('gray'); axis image; axis off;
title('bearbeitet');
