function WienerDeconvolution(NoiseLevel)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Demo der Wiener Entfaltung anhand von Platinen auf einem bewegten
%% Förderband.
%%
%% Das Eingangsbild x der unbewegten Platine wird Aufgrund der Bewegung des
%% Förderbandes während der Belichtung mit der zunächst unbekannten
%% Impulsantwort g gefaltet (Bewegungsunschärfe).
%% Zusätzlich sind die einzelnen Pixel mit Sensorrauschen behaftet, welches
%% für jedes Pixel als unabhängig von den anderen Pixeln angenommen wird.
%% Die Amplitude dieses Sensorrauschens sei normalverteilt.
%%
%% Wenn man eine Platine auf das gestoppte Förderband legt, und mehrere
%% Bilder aufnimmt, kann man über die Bilder mitteln um das Sensorrauschen
%% zu entfernen und so eine Schätzung des unverrauschten Bildes erhalten.
%%
%% Mit Hilfe dieses Bildes, und mehreren Bildern, welche bei 
%% eingeschaltetem Förderband aufgenommen wurden, lässt sich dann die
%% Impulsantwort g durch Systemidentifikation schätzen.
%%
%% Mit Hilfe der bekannten Impulsantwort g und den durch Periodogrammen
%% abschätzbaren Leistungsdichten von unverrauschtem Bild und Rauschen kann
%% über die Wiener Entfaltung das Rekonstruktionsfilter H entworfen werden,
%% welches bei laufendem Förderband aus den verschwommenen und verrauschten
%% Kamerabildern wieder ein unverrauschtes und scharfes Bild optimal 
%% schätzt.
%%
%% Eingabewerte:
%%      NoiseLevel kann die Werte 'none', 'low', 'medium', oder 'high' 
%%      haben. Default: 'low'
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inFileN = 30; %%Anzahl an Realisierungen des Rauschprozesses über die 
              %%gemittelt werden kann/soll
inFileExt = '.png';


%%Eingabewert überprüfen/auswerten
if nargin < 1 || isempty(NoiseLevel)
    NoiseLevel = 'low';
end
switch (NoiseLevel)
    case 'none'
        sig = 0;    
    case 'low'
        sig = 0.0001;
    case 'medium'
        sig = 0.001;
    case 'high'
        sig = 0.01;
    otherwise
        sig = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Verrauschte Eingangsdaten erzeugen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clean = imread('Platine_klein.png'); %Bild einlesen
name = 'Platine';

clean_moving = imread('Platine_klein_bewegt_clean.png'); %Bild einlesen
name_moving = 'Platine_bewegt';

for k=1:inFileN 
    noisy = imnoise(clean, 'gaussian', 0, sig); %Bild mit Rauschen versehen
    namek = strcat(name, num2str(k), '.png'); %Name des Ausgabe-Files erstellen
    imwrite(noisy, namek, 'png'); %Verrauschtes Bild speichern
    
    
    noisy = imnoise(clean_moving, 'gaussian', 0, sig); %Bild mit Rauschen versehen
    namek = strcat(name_moving, num2str(k), '.png'); %Name des Ausgabe-Files erstellen
    imwrite(noisy, namek', 'png'); %Verrauschtes Bild speichern
end
%%Variablen löschen um Arbeitsspeicher freizugeben
clear('noisy');
clear('namek');
clear('clean');
clear('clean_moving');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Sensorrauschen der Kamera schätzen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Es wird angenommen, dass ein Werkstück mehrmals in der gleichen Pose
%aufgenommen wurde. Da das Sensorrauschen als mittelwertfrei
%angenommen werden kann, betrachten wir das gemittelte Bild als den
%"wahren" Eingangswert x, und die Differenzen der einzelnen Bilder zum
%Mittelwertbild als Rauschen.

inFileName = 'Platine';
%%verrauschte Bilder ohne Bewegung öffnen
for k=1:inFileN
    filename = strcat(inFileName, num2str(k), inFileExt);
    %Speichere die Bilder in einem 4-D Bildstapel
    % 1. Dimension    : Bild-Nummer
    % 2.+3. Dimension : X/Y-Achse des Bildes
    % 4. Dimension    : Farbkanal (R=1, G=2, B=3)
    y1(k, :, :, :) = imread(filename); %Bezeichnung y1, da diese Bilder 
                                       %ohne Bewegung, d.h. für G = 1
                                       %aufgenommen wurden
end
[~, Y, X, ~] = size(y1); %%Bildgröße auslesen

%Mittelwert über alle Bilder bilden, also entlang der ersten
%Dimension, die damit überflüssig wird.
%squeeze() unterdrückt die unnötig gewordene Dimension, so dass
%x wieder ein Bild mit drei Kanälen ist.
x = squeeze(mean(y1, 1)); 

n = zeros(size(y1));
for k=1:inFileN
    n(k, :, :, :) = double(squeeze(y1(k,:,:,:))) - x; %Rauschen ist Differenzbild
end


%%Ausgabe der Zwischenergebnisse
figure;
imshow(uint8(x));
title('Geschätztes Eingangsbild ohne Rauschen');
figure;
imshow(uint8(squeeze(y1(2,:,:,:))));
title('Ein verrauschtes Eingangsbild');
figure;
imshow(uint8(squeeze(n(2,:,:,:))));
title('Geschätztes Rauschen');

%Haltepunkt zum Anschauen der Zwischenergebnisse
keyboard;
close all;

%%Variablen löschen um Arbeitsspeicher freizugeben
clear('y1');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%"Übertragunsfunktion" der Bewegung schätzen, Systemidentifikation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inFileName = 'Platine_bewegt';
%%verrauschte Bilder mit Bewegung öffnen
for k=1:inFileN
    filename = strcat(inFileName, num2str(k), inFileExt);
    %Speichere die Bilder in einem 4-D Bildstapel
    % 1. Dimension    : Bild nummer
    % 2.+3. Dimension : X/Y Achse des Bildes
    % 4. Dimension    : Farbkanal (R=1, G=2, B=3)
    y(k, :, :, :) = imread(filename);
end

%%Bilder enthalten noch Teile, die nicht zum Laufband gehören. Diese würden
%%unsere Schätzung verfälschen, darum betrachten wir nur einen vertikalen 
%%Streifen der Bilder (Zusatz _c für cropped)
x_c = x(floor(Y/4):ceil(3*Y/4), :, :);
y_c = y(:, floor(Y/4):ceil(3*Y/4), :, :);
n_c = n(:, floor(Y/4):ceil(3*Y/4), :, :);
[Y_c, X_c, ~] = size(x_c);

%%Leistungsdichtespektren von Ein- und Ausgangsignal schätzen. Wird über
%%Periodogramm gemacht
S_xx = zeros(Y_c, X_c);
for c=1:3 %%Mittelung über alle Farbkanäle
    S_xx = S_xx + 1/3*Periodogram2d( double(x_c(:,:,c)) );
end

S_yx = zeros(Y_c, X_c);
for k=1:inFileN %%Mittelung über alle Bilder
    for c=1:3 %%Mittelung über alle Farbkanäle
        S_yx = S_yx + 1/inFileN*1/3*Periodogram2d( double(squeeze(y_c(k,:,:,c))), x_c(:,:,c) );
    end
end

%Systemidentifikation über Leistungsdichteschätzungen
G = S_yx./S_xx;
g = fftshift(ifft2(G)); %Impulsantwort, in der Bildverarbeitung auch 
                        %point-spread function genannt
                        
%%Nur wenige Koeffizienten von g sind tatsächlich >> 0, so dass nur ein
%%Ausschnitt benötigt wird der jetzt ausgeschnitten wird
N_x = 20;
N_y = 3;
g_c = g( round(Y_c/2)-N_y:round(Y_c/2)+N_y, round(X_c/2)-N_x:round(X_c/2)+N_x);

%%Ausgabe der Zwischenergebnisse
figure;
imshow( uint8(x_c) );
title('Geschätztes Eingangsbild (ohne Bewegung), zugeschnitten');
figure;
imshow( squeeze(y_c(2,:,:,:)) );
title('Ausgangsbild (mit Bewegung und Rauschen), zugeschnitten');
figure;
imshow( uint8(double(intmax('uint8'))*g./(max(max(g)))) );
title('Geschätzte Impulsantwort der Bewegung');
figure;
imshow( uint8(double(intmax('uint8'))*g_c./(max(max(g_c)))) );
title('Geschätzte Impulsantwort, zugeschnitten');

%%Haltepunkt zum Anschauen der Zwischenergebnisse
keyboard;
close all;

%%Löschen von Variablen um Arbeitsspeicher freizugeben
clear('g');
clear('g_c');
clear('y');
clear('n');
clear('x');
%clear('S_yx');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Wiener Entfaltung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S_yy = zeros(Y_c, X_c);
%S_nn = zeros(Y_c, X_c);
for k=1:inFileN %%Mittelung über alle Bilder
    for c=1:3 %%Mittelung über alle Farbkanäle
        S_yy = S_yy + 1/inFileN*1/3*Periodogram2d( double(squeeze(y_c(k,:,:,c))) );
        %S_nn = S_nn + 1/inFileN*1/3*Periodogram2d( double(squeeze(n_c(k,:,:,c))) );
    end
end

%Hinweis:
%S_xy = conj(S_yx)

H = conj(S_yx)./S_yy;

%%alternative Berechnungsart, führt zum gleichen Ergebniss, allerdings muss 
%%S_nn anstatt S_yy explizit bestimmt werden (siehe ausgeklammerte Zeilen 
%%oben).
%H = S_xx.*conj(G)./(S_xx.*G.*conj(G)+S_nn);  

%%Faltung aller Eingangsbilder mit dem Rekonstruktionsfilter (aus
%%Geschwindigketsgründen im Frequenzbereich)
Y = zeros(size(y_c));
x_hat = zeros(size(y_c));
for k=1:inFileN %%Entfaltung für alle Bilder
    for c=1:3 %%Entfaltung für alle Farbkanäle, durch Faltung mit h
        Y(k, :, :, c) = fft2( double(squeeze(y_c(k,:,:,c))) ) ;
        x_hat(k, :, :, c) = ifft2( squeeze(Y(k, :, :, c)) .* H );
    end
end


%Ausgabe von (Zwischen-)Ergebnissen
figure
imshow( fftshift( uint8( double(intmax('uint8'))*abs( H )./max(max(abs(H))) ) ) );
title('|H(f)|');

figure
imshow( fftshift( uint8( double(intmax('uint8'))*abs( G )./max(max(abs(G))) ) ) );
title('|G(f)|');

figure
imshow( fftshift( uint8( double(intmax('uint8'))*abs(1/X_c/Y_c*squeeze(Y(2,:,:,1))) ) ) );
title('|Y(f)|');

figure
imshow( uint8(squeeze(y_c(2, :, :, :))) );
title('Eingangs Bild, zugeschnitten');

figure
imshow( uint8(squeeze(x_hat(2, :, :, :))) );
title('Rekonstruiertes Bild, zugeschnitten');

figure
imshow( uint8( x_c ) );
title('Korrektes Bild, zugeschnitten');

%Haltepunkt vor dem Ende der Demo
keyboard;
close all;


%Hilfsfunktion um die Leistungsdichte über Periodogramme zu schätzen
function LDS = Periodogram2d(A, B)
    fftA = fft2(A);
    if nargin < 2
        fftB = fftA;
    else
        fftB = fft2(B);
    end
    LDS = fftA .* conj(fftB);
end

end