% Luminanzkorrektur
%
% Aus gemessenen Luminanzwerten an verschiedenen Bildschirmpunkten mit
% verschiedenen Grauwerten werden für alle Punkte die Gamma-Funktionen
% (Abhängigkeit Luminanz voim Grauwert) quadratisch gefittet und eine
% Korrekturmatrix so errechnet, dass alle Luminanzen möglichst auf den
% dunkelsten Punkt angegelichen werden. Da die Korrektur konstant im Bezug
% auf den Grauwert ist, muss ein Luminanzwert angegeben werden, bei dem die
% Kompensation optimal ist.

clear all

% Konstanten
% gewünschte Bildschirmauflösung
X_Pixel = 1280;
Y_Pixel = 1024;
% Luminanzwert, bei dem die Kompensation optimal ist
SollLuminanz = 75; % cd/m2

% Datei mit Kalibrierungsdaten
filename='merged_calibrate_display.2008-02-19_13_48_52_calibrate_display.2008-02-21_08_57_08.cal';

% Daten lesen
[ mode	pos	x	y	grey	lum ] = textread( filename, '%d%d%d%d%f%f', 'headerlines', 13);

% Verteilung bei Grauwert 128
% Zeilen (Y-Werte)
N=15;
% Spalten (X-Werte)
M=20;
% hier kommen die gemessenen Luminanzwerte rein (2 Messungen)
lumGes1 = zeros(N,M);
lumGes2 = zeros(N,M);
% Modus 1: Messung mit hoher Auflösung
modeVal = 1;
idx = find(mode==modeVal);
% in x und y stehen die Positionen an denen gemessen wurde.
% (0,0)=links unten
x128=x(idx);
y128=y(idx);
% grober plot
%plot3(x(idx),y(idx),lum(idx),'.')
% in due beiden lumGes-Matrizen einsortieren
for i=1:length(idx)
    index=idx(i);
    position=pos(index);
    n=rem(position,M)+1;
    m=floor(position/M)+1;
    % erste Messung
    if lumGes1(m,n)==0
        lumGes1(m,n)=lum(index);
    % zweite Messung
    else
        lumGes2(m,n)=lum(index);
    end
end
% darstellen
figure(1)
clf
subplot(2,2,1)
pcolor(lumGes1)
shading flat
title 'Luminanzen 128 erste Messung'
subplot(2,2,3)
surf(lumGes1)
subplot(2,2,2)
pcolor(lumGes2)
shading flat
title 'Luminanzen 128 zweite Messung'
subplot(2,2,4)
surf(lumGes2)
% Differenz der Messungen darstellen
figure(2)
surf(lumGes1-lumGes2)
l1=lumGes1(:);
l2=lumGes2(:);
mittlAbw = mean(l2-l1);
mittlAbwStd = std(l2-l1);
title(['mittlere Abweichung: ' num2str(mittlAbw) '+-' num2str(mittlAbwStd) ' cd/m2'])
% Mittelwert
lumGes = (lumGes1+lumGes2)/2;
lumGes128 = lumGes;

% Gamma-Kurven
% Modus ist hier 0
modeVal = 0;
idx = find(mode==modeVal);
% Positionen
x=x(idx);
y=y(idx);
pos=pos(idx);
% Grauwerte
grey=grey(idx);
% Luminanzen
lum=lum(idx);
% Gammafunktionen (Luminanzen je Grauwert und Position)
gamma=zeros(10,9);
% Positionen (9 Stück = 3x3)
xgrey=zeros(9,1);
ygrey=zeros(9,1);
figure(3)
clf
% Es sind je Position 10 Grauwerte gemessen worden
for greyValIdx=1:10
    N=3;
    M=3;
    % für jede der 2 Messungen eine Luminanzmatrix
    lumGes1 = ones(N,M).*-100;
    lumGes2 = ones(N,M).*-100;
    % die Grauwerte gehen von 0 bis 255
    greyVal(greyValIdx) =(greyValIdx-1)*(255/9);
    % entsprechende Werte aus den Messdaten herausfischen
    idxGreyVal = find(abs(grey-greyVal(greyValIdx))<0.01);
    % für alle Grauwerte-Messungen
    for i=1:length(idxGreyVal)
        index=idxGreyVal(i);
        % Positionsindex
        position=pos(index);
        % x-y-Positionen
        xgrey(position+1) = x(index);
        ygrey(position+1) = y(index);
        % Position in der Matrix berechnen
        n=rem(position,M)+1;
        m=floor(position/M)+1;
        % erste Messung
        if lumGes1(m,n)==-100
            lumGes1(m,n)=lum(index);
            gamma(greyValIdx, position+1)=lum(index);
        else
        % zweite Messung, Luminanz mitteln
            lumGes2(m,n)=lum(index);
            gamma(greyValIdx, position+1)= (gamma(greyValIdx, position+1)+ lum(index))/2;
        end
    end
    % darstellen
    subplot(2,2,1)
    hold on
    surf(lumGes1)
    subplot(2,2,2)
    hold on
    surf(lumGes2)
end

% Gamma-Funktionen plotten
figure(4)
cla
hold on
plot(gamma,'r')

% Korrektur für die 3x3-Daten aus Modus 0 (hochauflösende Korrektur folgt
% unten):
% "Sonnenbrille" simulieren:
% berechne für alle Positionen den Grauwert, der für eine bestimmte Luminanz erforderlich ist
Luminanz = gamma;
Grauwert = greyVal';
fig=figure(5);
for Position=1:9
    coefficient4gamma(Position, :) = polyfit(Grauwert, Luminanz(:,Position), 2); % Calculate 2nd order fit
    clf;
    plot(Grauwert, Luminanz(:,Position), 'k')
    hold on
    xnew = linspace(min(Grauwert), max(Grauwert), 100);
    plot(xnew, polyval(coefficient4gamma(Position,:), xnew), 'r')
    a=axis; axis([min(Grauwert) max(Grauwert) a(3) a(4)]); clear a
    title(['Position: ' num2str(Position)])
    xlabel('Grey Value')
    ylabel('Luminance [cd/m^2]')
    drawnow
    hold off
    % Mitternachtsformel zum invertieren des Polynoms 2. Ordnung:
    % 0 = c(3) + x*c(2) + x^2*c(1) - L
    c=coefficient4gamma(Position, :);
    SollGrauwert(Position) = ( -c(2)+sqrt( c(2)^2 - 4*c(1)*(c(3)-SollLuminanz) ) ) / ( 2 * c(1) );
end % if, debug
close(fig)

% Faktoren für die Abdunklung auf den kleinsten Wert berechnen
% Alle Grauwerte sollen auf den höchsten angegelichen werden
Maxi = max(SollGrauwert);
FaktorenGrauwert = SollGrauwert./Maxi;
FaktorMatrixGrauwert = [ FaktorenGrauwert(1:3); 
                        FaktorenGrauwert(4:6); 
                        FaktorenGrauwert(7:9) ]; 

% Anwendung dieser Korrekturfaktoren simulieren:
gamma_korr=zeros(10,9);
figure(3)
subplot(2,2,3)
cla
hold on
KorrigierterGrauwert=zeros(10,9);
for greyValIdx=1:10
    % Simulation
    N=3;
    M=3;
    lumGes = ones(N,M).*-100;
    % tatsächliche Luminanzen berechnen
    for Position = 1:9
        % Luminanz für korrigierten Grauwert berechnen (Simulation)
        KorrigierterGrauwert(greyValIdx, Position) = Grauwert(greyValIdx)*FaktorenGrauwert(Position);
        IstLuminanz = polyval(coefficient4gamma(Position,:), KorrigierterGrauwert(greyValIdx, Position) );
        % Luminanz-Matrix bauen
        n=rem(Position-1,M)+1;
        m=floor((Position-1)/M)+1;
        % erste Messung
        lumGes(m,n)=IstLuminanz;   
        gamma_korr(greyValIdx,Position) = IstLuminanz;
    end
    surf(lumGes)
end
% resultierende korrigierte Gamma-Funktionen plotten
figure(4)
plot(gamma_korr,'k')

% Korrekturmatrix für alle Lum128-Werte erstellen.
% Dafür wird für jeden Punkt die gefittete mittlere gamma-Funktion entsprechend skaliert.
% mittlere gamma-Funktion
gamma_mean = mean(gamma');
% Calculate 3rd order fit
coefficient4gammamean = polyfit(Grauwert, gamma_mean', 2); 
%lum_mean_128 = polyval(coefficient4gammamean, 128 );

% Verteilung bei Grauwert 128
% Daten nochmal einlesen
[ mode	pos	x	y	grey	lum ] = textread( filename, '%d%d%d%d%f%f', 'headerlines', 13);
% Dimension der höher aufgelösten Messmatrix
N=15;
M=20;
% Matrizen für die Luminanzen initialisieren
lumGes1 = zeros(N,M);
lumGes2 = zeros(N,M);
% jetzt holen wir wieder die Daten vom Modus 1
modeVal = 1;
idx = find(mode==modeVal);
% für jeden Messwert
for i=1:length(idx)
    index=idx(i);
    position=pos(index);
    % Indizes der Matrix berechnen
    n=rem(position,M)+1;
    m=floor(position/M)+1;
    % erste Messung
    if lumGes1(m,n)==0
        lumGes1(m,n)=lum(index);
        % es reicht, die Positionen bei der ersten Messung zu holen
        x128(i)=x(index);
        y128(i)=y(index);
    % zweite Messung
    else
        lumGes2(m,n)=lum(index);
    end
end
% Luminanzen 1-dimensional machen und die beiden Messungen mitteln
Lum1 = lumGes1;
Lum1 = Lum1';
Lum1 = Lum1(:);
Lum2 = lumGes2;
Lum2 = Lum2';
Lum2 = Lum2(:);
Lum = mean([Lum1 Lum2]');
Lum = Lum';
% für alle Positionen
N=size(Lum,1);
for i=1:N
    % Verhältnis zum mittleren lum128-Wert bestimmen
    Faktor = Lum(i) / polyval(coefficient4gammamean, 128 );
    % gamma-Funktion entsprechend skalieren 
    % (die Koeffizienten, weil's ein Polynom ist)
    Koeffizienten(i,:)=coefficient4gammamean.*Faktor;
    % Soll-Grauwert (Sonnenbrille) berechnen
    c=Koeffizienten(i, :);
    SollGrauwert(i) = ( -c(2)+sqrt( c(2)^2 - 4*c(1)*(c(3)-SollLuminanz) ) ) / ( 2 * c(1) );
end
% Faktoren für die Abdunklung auf den kleinsten Wert berechnen
% Alle Grauwerte sollen auf den höchsten angegelichen werden
Maxi = max(SollGrauwert);
FaktorenGrauwert = SollGrauwert./Maxi;

% Simulation der Korrektur
KorrigierterGrauwert=zeros(N,1);
N=15;
M=20;
% Sollenbrillenmatrix
FaktorenGrauwertMatrix=zeros(N,M);
% Simulierte Luminanz-Matrix
lumGes = zeros(N,M);
% tatsächliche Luminanzen berechnen
for Position = 1:size(Lum,1)
    % Luminanz für korrigierten Grauwert berechnen (Simulation)
    KorrigierterGrauwert(Position) = 128*FaktorenGrauwert(Position);
    IstLuminanz = polyval(Koeffizienten(Position,:), KorrigierterGrauwert(Position) );
    % Luminanz-Matrix bauen und Index berechnen
    n=rem(Position-1,M)+1;
    m=floor((Position-1)/M)+1;
    % simulierten Luminanzwert abspeichern
    lumGes(m,n)=IstLuminanz;
    % Sollenbrillenmatrix: Wert eintragen
    FaktorenGrauwertMatrix(m,n)=FaktorenGrauwert(Position);
end
% Darstellung
figure(5)
clf
surf(lumGes)
hold on
surf(lumGes128,'FaceAlpha',0.2)

% mit interp2 volle Bildschirm-Pixel-KorrekturMatrix erzeugen
% Messungs-Positionen
xv=unique(x)';
yvals=unique(y)';
% verschiedene Positionen der beiden Messungen zusammenführen
xvals=xv;
xvals(11) = mean(xvals(11:12));
xvals(12)=[];
[xg,yg] = meshgrid(xvals,yvals);
% Pixelpositionen
[xi,yi] = meshgrid(1:X_Pixel,1:Y_Pixel);
%[xi,yi] = meshgrid(min(xvals):max(xvals),min(yvals):max(yvals));
% Interpolieren
% Ausserhalb des Messbereichs werden die Werte auf 0 gesetzt, weil sie
% sonst bei der Interpolation abhauen.
zi = interp2(xg,yg,FaktorenGrauwertMatrix,xi,yi,'linear',0); 
clf
% Plot-Daten für die Darstellung nochmal interpolieren (sind sonst zu viele)
[xn,yn] = meshgrid(1:20:X_Pixel,1:20:Y_Pixel);
zn=interp2(xi,yi,zi,xn,yn,'linear',0);
mesh(xn,yn,zn)
hold on
surf(xg,yg,FaktorenGrauwertMatrix,'FaceAlpha',0.2)
% in zi stehen die interpolierten Korrekturfaktoren
% Abspeichern in eine Datei, zeilenweise statt spaltenweise: Transponieren
% Linux-Rechner will ieee-le - Format
% GL_FLOAT für Alpha-Textur braucht 32-Bit floats
f=fopen('LuminanceCorrection1280x1024.dat','w','ieee-le');
fwrite(f,zi','float32');
fclose(f);

% zum Testen, ob die Orientierung stimmt:
a=zi;
for i=100:300
    for j= 400:700
        a(i,j)=0;
    end
end
figure(6)
zn=interp2(xi,yi,a,xn,yn,'linear',0);
mesh(xn,yn,zn)
f=fopen('LuminanceCorrection1280x1024Test.dat','w','ieee-le');
fwrite(f,a','float32');
fclose(f);

% Ergebnisse der Korrekturmessung
figure(8)
clf
surf(xg,yg,lumGes128)
hold on
grid on
for idx=2:3
    fname=sprintf('12x9_128_%d.txt',idx);
    lum_korr_gemessen=textread(fname);

    N=9;
    M=12;
    % Gemessene Luminanz-Matrix
    lumKorr = zeros(N,M);
    % tatsächliche Luminanzen berechnen
    for Position = 1:size(lum_korr_gemessen,1)
        % Luminanz-Matrix bauen und Index berechnen
        n=rem(Position-1,M)+1;
        m=floor((Position-1)/M)+1;
        % simulierten Luminanzwert abspeichern
        lumKorr(m,n)=lum_korr_gemessen( Position );
    end
    % Positionen in Pixel umrechnen
    % Bildschirm: 375x300 mm 1280x1024 Pixel Abstand 500 mm
    % in Grad
    xpos = linspace(-13.44,13.44,12)';
    ypos = linspace(-9.31,9.31,9)';
    xvals = 500 * tan( pi/180*xpos ) * 1280 / 375 + X_Pixel/2;
    yvals = 500 * tan( pi/180*ypos ) * 1024 / 300 + Y_Pixel/2;
    % Pixel-Matrix bauen
    [xg,yg] = meshgrid(xvals,yvals);
    % % Pixelpositionen
    % [xi,yi] = meshgrid(1:X_Pixel,1:Y_Pixel);
    % % Interpolieren
    % % Ausserhalb des Messbereichs werden die Werte auf 0 gesetzt, weil sie
    % % sonst bei der Interpolation abhauen.
    % zi = interp2(xg,yg,lumKorr,xi,yi,'linear',0);
    %
    % % Plot-Daten für die Darstellung nochmal interpolieren (sind sonst zu viele)
    % [xn,yn] = meshgrid(1:20:X_Pixel,1:20:Y_Pixel);
    % zn=interp2(xi,yi,zi,xn,yn,'linear',0);
    surf(xg,yg,lumKorr)
end
