1 % Luminanzkorrektur
   2 %
   3 % Aus gemessenen Luminanzwerten an verschiedenen Bildschirmpunkten mit
   4 % verschiedenen Grauwerten werden für alle Punkte die Gamma-Funktionen
   5 % (Abhängigkeit Luminanz voim Grauwert) quadratisch gefittet und eine
   6 % Korrekturmatrix so errechnet, dass alle Luminanzen möglichst auf den
   7 % dunkelsten Punkt angegelichen werden. Da die Korrektur konstant im Bezug
   8 % auf den Grauwert ist, muss ein Luminanzwert angegeben werden, bei dem die
   9 % Kompensation optimal ist.
  10 
  11 clear all
  12 
  13 % Konstanten
  14 % gewünschte Bildschirmauflösung
  15 X_Pixel = 1280;
  16 Y_Pixel = 1024;
  17 % Luminanzwert, bei dem die Kompensation optimal ist
  18 SollLuminanz = 75; % cd/m2
  19 
  20 % Datei mit Kalibrierungsdaten
  21 filename='merged_calibrate_display.2008-02-19_13_48_52_calibrate_display.2008-02-21_08_57_08.cal';
  22 
  23 % Daten lesen
  24 [ mode	pos	x	y	grey	lum ] = textread( filename, '%d%d%d%d%f%f', 'headerlines', 13);
  25 
  26 % Verteilung bei Grauwert 128
  27 % Zeilen (Y-Werte)
  28 N=15;
  29 % Spalten (X-Werte)
  30 M=20;
  31 % hier kommen die gemessenen Luminanzwerte rein (2 Messungen)
  32 lumGes1 = zeros(N,M);
  33 lumGes2 = zeros(N,M);
  34 % Modus 1: Messung mit hoher Auflösung
  35 modeVal = 1;
  36 idx = find(mode==modeVal);
  37 % in x und y stehen die Positionen an denen gemessen wurde.
  38 % (0,0)=links unten
  39 x128=x(idx);
  40 y128=y(idx);
  41 % grober plot
  42 %plot3(x(idx),y(idx),lum(idx),'.')
  43 % in due beiden lumGes-Matrizen einsortieren
  44 for i=1:length(idx)
  45     index=idx(i);
  46     position=pos(index);
  47     n=rem(position,M)+1;
  48     m=floor(position/M)+1;
  49     % erste Messung
  50     if lumGes1(m,n)==0
  51         lumGes1(m,n)=lum(index);
  52     % zweite Messung
  53     else
  54         lumGes2(m,n)=lum(index);
  55     end
  56 end
  57 % darstellen
  58 figure(1)
  59 clf
  60 subplot(2,2,1)
  61 pcolor(lumGes1)
  62 shading flat
  63 title 'Luminanzen 128 erste Messung'
  64 subplot(2,2,3)
  65 surf(lumGes1)
  66 subplot(2,2,2)
  67 pcolor(lumGes2)
  68 shading flat
  69 title 'Luminanzen 128 zweite Messung'
  70 subplot(2,2,4)
  71 surf(lumGes2)
  72 % Differenz der Messungen darstellen
  73 figure(2)
  74 surf(lumGes1-lumGes2)
  75 l1=lumGes1(:);
  76 l2=lumGes2(:);
  77 mittlAbw = mean(l2-l1);
  78 mittlAbwStd = std(l2-l1);
  79 title(['mittlere Abweichung: ' num2str(mittlAbw) '+-' num2str(mittlAbwStd) ' cd/m2'])
  80 % Mittelwert
  81 lumGes = (lumGes1+lumGes2)/2;
  82 lumGes128 = lumGes;
  83 
  84 % Gamma-Kurven
  85 % Modus ist hier 0
  86 modeVal = 0;
  87 idx = find(mode==modeVal);
  88 % Positionen
  89 x=x(idx);
  90 y=y(idx);
  91 pos=pos(idx);
  92 % Grauwerte
  93 grey=grey(idx);
  94 % Luminanzen
  95 lum=lum(idx);
  96 % Gammafunktionen (Luminanzen je Grauwert und Position)
  97 gamma=zeros(10,9);
  98 % Positionen (9 Stück = 3x3)
  99 xgrey=zeros(9,1);
 100 ygrey=zeros(9,1);
 101 figure(3)
 102 clf
 103 % Es sind je Position 10 Grauwerte gemessen worden
 104 for greyValIdx=1:10
 105     N=3;
 106     M=3;
 107     % für jede der 2 Messungen eine Luminanzmatrix
 108     lumGes1 = ones(N,M).*-100;
 109     lumGes2 = ones(N,M).*-100;
 110     % die Grauwerte gehen von 0 bis 255
 111     greyVal(greyValIdx) =(greyValIdx-1)*(255/9);
 112     % entsprechende Werte aus den Messdaten herausfischen
 113     idxGreyVal = find(abs(grey-greyVal(greyValIdx))<0.01);
 114     % für alle Grauwerte-Messungen
 115     for i=1:length(idxGreyVal)
 116         index=idxGreyVal(i);
 117         % Positionsindex
 118         position=pos(index);
 119         % x-y-Positionen
 120         xgrey(position+1) = x(index);
 121         ygrey(position+1) = y(index);
 122         % Position in der Matrix berechnen
 123         n=rem(position,M)+1;
 124         m=floor(position/M)+1;
 125         % erste Messung
 126         if lumGes1(m,n)==-100
 127             lumGes1(m,n)=lum(index);
 128             gamma(greyValIdx, position+1)=lum(index);
 129         else
 130         % zweite Messung, Luminanz mitteln
 131             lumGes2(m,n)=lum(index);
 132             gamma(greyValIdx, position+1)= (gamma(greyValIdx, position+1)+ lum(index))/2;
 133         end
 134     end
 135     % darstellen
 136     subplot(2,2,1)
 137     hold on
 138     surf(lumGes1)
 139     subplot(2,2,2)
 140     hold on
 141     surf(lumGes2)
 142 end
 143 
 144 % Gamma-Funktionen plotten
 145 figure(4)
 146 cla
 147 hold on
 148 plot(gamma,'r')
 149 
 150 % Korrektur für die 3x3-Daten aus Modus 0 (hochauflösende Korrektur folgt
 151 % unten):
 152 % "Sonnenbrille" simulieren:
 153 % berechne für alle Positionen den Grauwert, der für eine bestimmte Luminanz erforderlich ist
 154 Luminanz = gamma;
 155 Grauwert = greyVal';
 156 fig=figure(5);
 157 for Position=1:9
 158     coefficient4gamma(Position, :) = polyfit(Grauwert, Luminanz(:,Position), 2); % Calculate 2nd order fit
 159     clf;
 160     plot(Grauwert, Luminanz(:,Position), 'k')
 161     hold on
 162     xnew = linspace(min(Grauwert), max(Grauwert), 100);
 163     plot(xnew, polyval(coefficient4gamma(Position,:), xnew), 'r')
 164     a=axis; axis([min(Grauwert) max(Grauwert) a(3) a(4)]); clear a
 165     title(['Position: ' num2str(Position)])
 166     xlabel('Grey Value')
 167     ylabel('Luminance [cd/m^2]')
 168     drawnow
 169     hold off
 170     % Mitternachtsformel zum invertieren des Polynoms 2. Ordnung:
 171     % 0 = c(3) + x*c(2) + x^2*c(1) - L
 172     c=coefficient4gamma(Position, :);
 173     SollGrauwert(Position) = ( -c(2)+sqrt( c(2)^2 - 4*c(1)*(c(3)-SollLuminanz) ) ) / ( 2 * c(1) );
 174 end % if, debug
 175 close(fig)
 176 
 177 % Faktoren für die Abdunklung auf den kleinsten Wert berechnen
 178 % Alle Grauwerte sollen auf den höchsten angegelichen werden
 179 Maxi = max(SollGrauwert);
 180 FaktorenGrauwert = SollGrauwert./Maxi;
 181 FaktorMatrixGrauwert = [ FaktorenGrauwert(1:3); 
 182                         FaktorenGrauwert(4:6); 
 183                         FaktorenGrauwert(7:9) ]; 
 184 
 185 % Anwendung dieser Korrekturfaktoren simulieren:
 186 gamma_korr=zeros(10,9);
 187 figure(3)
 188 subplot(2,2,3)
 189 cla
 190 hold on
 191 KorrigierterGrauwert=zeros(10,9);
 192 for greyValIdx=1:10
 193     % Simulation
 194     N=3;
 195     M=3;
 196     lumGes = ones(N,M).*-100;
 197     % tatsächliche Luminanzen berechnen
 198     for Position = 1:9
 199         % Luminanz für korrigierten Grauwert berechnen (Simulation)
 200         KorrigierterGrauwert(greyValIdx, Position) = Grauwert(greyValIdx)*FaktorenGrauwert(Position);
 201         IstLuminanz = polyval(coefficient4gamma(Position,:), KorrigierterGrauwert(greyValIdx, Position) );
 202         % Luminanz-Matrix bauen
 203         n=rem(Position-1,M)+1;
 204         m=floor((Position-1)/M)+1;
 205         % erste Messung
 206         lumGes(m,n)=IstLuminanz;   
 207         gamma_korr(greyValIdx,Position) = IstLuminanz;
 208     end
 209     surf(lumGes)
 210 end
 211 % resultierende korrigierte Gamma-Funktionen plotten
 212 figure(4)
 213 plot(gamma_korr,'k')
 214 
 215 % Korrekturmatrix für alle Lum128-Werte erstellen.
 216 % Dafür wird für jeden Punkt die gefittete mittlere gamma-Funktion entsprechend skaliert.
 217 % mittlere gamma-Funktion
 218 gamma_mean = mean(gamma');
 219 % Calculate 3rd order fit
 220 coefficient4gammamean = polyfit(Grauwert, gamma_mean', 2); 
 221 %lum_mean_128 = polyval(coefficient4gammamean, 128 );
 222 
 223 % Verteilung bei Grauwert 128
 224 % Daten nochmal einlesen
 225 [ mode	pos	x	y	grey	lum ] = textread( filename, '%d%d%d%d%f%f', 'headerlines', 13);
 226 % Dimension der höher aufgelösten Messmatrix
 227 N=15;
 228 M=20;
 229 % Matrizen für die Luminanzen initialisieren
 230 lumGes1 = zeros(N,M);
 231 lumGes2 = zeros(N,M);
 232 % jetzt holen wir wieder die Daten vom Modus 1
 233 modeVal = 1;
 234 idx = find(mode==modeVal);
 235 % für jeden Messwert
 236 for i=1:length(idx)
 237     index=idx(i);
 238     position=pos(index);
 239     % Indizes der Matrix berechnen
 240     n=rem(position,M)+1;
 241     m=floor(position/M)+1;
 242     % erste Messung
 243     if lumGes1(m,n)==0
 244         lumGes1(m,n)=lum(index);
 245         % es reicht, die Positionen bei der ersten Messung zu holen
 246         x128(i)=x(index);
 247         y128(i)=y(index);
 248     % zweite Messung
 249     else
 250         lumGes2(m,n)=lum(index);
 251     end
 252 end
 253 % Luminanzen 1-dimensional machen und die beiden Messungen mitteln
 254 Lum1 = lumGes1;
 255 Lum1 = Lum1';
 256 Lum1 = Lum1(:);
 257 Lum2 = lumGes2;
 258 Lum2 = Lum2';
 259 Lum2 = Lum2(:);
 260 Lum = mean([Lum1 Lum2]');
 261 Lum = Lum';
 262 % für alle Positionen
 263 N=size(Lum,1);
 264 for i=1:N
 265     % Verhältnis zum mittleren lum128-Wert bestimmen
 266     Faktor = Lum(i) / polyval(coefficient4gammamean, 128 );
 267     % gamma-Funktion entsprechend skalieren 
 268     % (die Koeffizienten, weil's ein Polynom ist)
 269     Koeffizienten(i,:)=coefficient4gammamean.*Faktor;
 270     % Soll-Grauwert (Sonnenbrille) berechnen
 271     c=Koeffizienten(i, :);
 272     SollGrauwert(i) = ( -c(2)+sqrt( c(2)^2 - 4*c(1)*(c(3)-SollLuminanz) ) ) / ( 2 * c(1) );
 273 end
 274 % Faktoren für die Abdunklung auf den kleinsten Wert berechnen
 275 % Alle Grauwerte sollen auf den höchsten angegelichen werden
 276 Maxi = max(SollGrauwert);
 277 FaktorenGrauwert = SollGrauwert./Maxi;
 278 
 279 % Simulation der Korrektur
 280 KorrigierterGrauwert=zeros(N,1);
 281 N=15;
 282 M=20;
 283 % Sollenbrillenmatrix
 284 FaktorenGrauwertMatrix=zeros(N,M);
 285 % Simulierte Luminanz-Matrix
 286 lumGes = zeros(N,M);
 287 % tatsächliche Luminanzen berechnen
 288 for Position = 1:size(Lum,1)
 289     % Luminanz für korrigierten Grauwert berechnen (Simulation)
 290     KorrigierterGrauwert(Position) = 128*FaktorenGrauwert(Position);
 291     IstLuminanz = polyval(Koeffizienten(Position,:), KorrigierterGrauwert(Position) );
 292     % Luminanz-Matrix bauen und Index berechnen
 293     n=rem(Position-1,M)+1;
 294     m=floor((Position-1)/M)+1;
 295     % simulierten Luminanzwert abspeichern
 296     lumGes(m,n)=IstLuminanz;
 297     % Sollenbrillenmatrix: Wert eintragen
 298     FaktorenGrauwertMatrix(m,n)=FaktorenGrauwert(Position);
 299 end
 300 % Darstellung
 301 figure(5)
 302 clf
 303 surf(lumGes)
 304 hold on
 305 surf(lumGes128,'FaceAlpha',0.2)
 306 
 307 % mit interp2 volle Bildschirm-Pixel-KorrekturMatrix erzeugen
 308 % Messungs-Positionen
 309 xv=unique(x)';
 310 yvals=unique(y)';
 311 % verschiedene Positionen der beiden Messungen zusammenführen
 312 xvals=xv;
 313 xvals(11) = mean(xvals(11:12));
 314 xvals(12)=[];
 315 [xg,yg] = meshgrid(xvals,yvals);
 316 % Pixelpositionen
 317 [xi,yi] = meshgrid(1:X_Pixel,1:Y_Pixel);
 318 %[xi,yi] = meshgrid(min(xvals):max(xvals),min(yvals):max(yvals));
 319 % Interpolieren
 320 % Ausserhalb des Messbereichs werden die Werte auf 0 gesetzt, weil sie
 321 % sonst bei der Interpolation abhauen.
 322 zi = interp2(xg,yg,FaktorenGrauwertMatrix,xi,yi,'linear',0); 
 323 clf
 324 % Plot-Daten für die Darstellung nochmal interpolieren (sind sonst zu viele)
 325 [xn,yn] = meshgrid(1:20:X_Pixel,1:20:Y_Pixel);
 326 zn=interp2(xi,yi,zi,xn,yn,'linear',0);
 327 mesh(xn,yn,zn)
 328 hold on
 329 surf(xg,yg,FaktorenGrauwertMatrix,'FaceAlpha',0.2)
 330 % in zi stehen die interpolierten Korrekturfaktoren
 331 % Abspeichern in eine Datei, zeilenweise statt spaltenweise: Transponieren
 332 % Linux-Rechner will ieee-le - Format
 333 % GL_FLOAT für Alpha-Textur braucht 32-Bit floats
 334 f=fopen('LuminanceCorrection1280x1024.dat','w','ieee-le');
 335 fwrite(f,zi','float32');
 336 fclose(f);
 337 
 338 % zum Testen, ob die Orientierung stimmt:
 339 a=zi;
 340 for i=100:300
 341     for j= 400:700
 342         a(i,j)=0;
 343     end
 344 end
 345 figure(6)
 346 zn=interp2(xi,yi,a,xn,yn,'linear',0);
 347 mesh(xn,yn,zn)
 348 f=fopen('LuminanceCorrection1280x1024Test.dat','w','ieee-le');
 349 fwrite(f,a','float32');
 350 fclose(f);
 351 
 352 % Ergebnisse der Korrekturmessung
 353 figure(8)
 354 clf
 355 surf(xg,yg,lumGes128)
 356 hold on
 357 grid on
 358 for idx=2:3
 359     fname=sprintf('12x9_128_%d.txt',idx);
 360     lum_korr_gemessen=textread(fname);
 361 
 362     N=9;
 363     M=12;
 364     % Gemessene Luminanz-Matrix
 365     lumKorr = zeros(N,M);
 366     % tatsächliche Luminanzen berechnen
 367     for Position = 1:size(lum_korr_gemessen,1)
 368         % Luminanz-Matrix bauen und Index berechnen
 369         n=rem(Position-1,M)+1;
 370         m=floor((Position-1)/M)+1;
 371         % simulierten Luminanzwert abspeichern
 372         lumKorr(m,n)=lum_korr_gemessen( Position );
 373     end
 374     % Positionen in Pixel umrechnen
 375     % Bildschirm: 375x300 mm 1280x1024 Pixel Abstand 500 mm
 376     % in Grad
 377     xpos = linspace(-13.44,13.44,12)';
 378     ypos = linspace(-9.31,9.31,9)';
 379     xvals = 500 * tan( pi/180*xpos ) * 1280 / 375 + X_Pixel/2;
 380     yvals = 500 * tan( pi/180*ypos ) * 1024 / 300 + Y_Pixel/2;
 381     % Pixel-Matrix bauen
 382     [xg,yg] = meshgrid(xvals,yvals);
 383     % % Pixelpositionen
 384     % [xi,yi] = meshgrid(1:X_Pixel,1:Y_Pixel);
 385     % % Interpolieren
 386     % % Ausserhalb des Messbereichs werden die Werte auf 0 gesetzt, weil sie
 387     % % sonst bei der Interpolation abhauen.
 388     % zi = interp2(xg,yg,lumKorr,xi,yi,'linear',0);
 389     %
 390     % % Plot-Daten für die Darstellung nochmal interpolieren (sind sonst zu viele)
 391     % [xn,yn] = meshgrid(1:20:X_Pixel,1:20:Y_Pixel);
 392     % zn=interp2(xi,yi,zi,xn,yn,'linear',0);
 393     surf(xg,yg,lumKorr)
 394 end