Úloha è

Úloha č.2
Zadání

Proveďte projektivní rekonstrukci ze dvou zadaných obrázků pořízených jednou kamerou.

Výpočet epipolární geometrie

Ze souřadnic bodů vyplníme matici A o rozměrech 9xn kde n je počet bodů:

A.q=0

q je sloupcový vektor který obsahuje prvky matice Q.

Matici Q získáme jako nulový prostor matice A.

Výpočet koeficientů lambda vychází ze vztahu:

Z koeficientů lambda snadno zobrazíme epipoláry:

lambda = {a, b, c} : a*x + b*y + c = 0

Vypočet epipólů:

Epipóly získáme z matice Q pomocí SVD, jako sloupce matic U a V.

Obr 1: Znázorňuje epipoláry v prvním obrázku Obr 2: Znázorňuje epipoláry v druhém obrázku

Původní zadané body jsou zobrazeny modře.

Obr 3: Obrázek znázorňuje epipól (červeně) v prvním obrázku. V druhém obrázku vyšel epipól tak daleko, že se nám ho nepodařilo zobrazit.

Rekonstrukce z kalibrovaných kamer

S využitím známých parametrů kamer ze souboru pointcal_mat rekonstruujeme 3D body z rovnic:

Obr 4: Rekonstrukce 3D Obr 4: Rekonstrukce 3D – jiný pohled

Rekonstruované body pomocí známých parametrů kamer se jsou znázorněny modře.

Z rekonstruovaných bodů provedeme reprojekci do obou obrazů. Ta je zobrazena na Obr 1 a 2 žlutou barvou.

Projektivni rekonstrukce

Zvolíme matici myP1=[O3 I3]. K ní vypočítáme matici myP2 pomocí matice W. Parametry matice W získáme z matice Q.

Matice W má hodnost 5.

Nalezneme bázi matice W. Ta má rozměr 9x4. Potom volíme lineární kombinaci vektorů báze matice W takovou, která je nejméně singulární. Tu najdeme pomocí poměru prvního a posledního čísla s , který musí být co největší.

Matici a potom získáme přerovnáním prvků výše zmíněného vektoru do matice 3x3. Z této matice vytvoříme matici :

myP2=[ A b ]

Po té rekonstruujeme 3D body pomocí matic myP1 a myP2. Provedeme reprojekci rekonstruovaných bodů do původních obrázků podle výše zmíněné postupu. Body jsou zobrazeny na Obr 1 a Obr 2 tyrkysovou barvou.

Takto rekonstruované 3D body se liší od skutečných bodů projektivní transformací. Abychom tuto transformaci odstranili použijeme matici H.

matici H a konstanty m i získáme jako nulový prostor matice Z, která má rozměr 4*n x 5*n kde n je počet bodů.

Nulový prostor má dimenzi 21, kde prvních n prvků jsou konstanty m i a posledních 3*n prvků jsou prvky matice H.

Pomocí matice H a známých konstant m i spočítáme rekonstruované 3D body.

Rekonstruované 3D body jsou zobrazeny v Obr 4 a Obr 5 zeleně. Dále provedeme reprojekci rekonstruovaných bodů do původních obrázků podle výše zmíněné postupu. Body jsou zobrzeny na Obr 1 a Obr 2 fialovou barvou.

Závěr

Existence rekonstrukce je podmíněna tím, že zadané body neleźí v rovině. Čím více je rozložení bodů rovinné, tím citlivější je rekonstukce na změnu polohy jednotlivých korespondujících bodů a poměrně malá chyba při označování korespondujících bodů může způsobit velkou chybu při výpočtu matice H. Mezi námi zvolenými body však lze proložit rovinu, od které nemají body příliš velikou odchylku. To pravděpodobně způsobilo větší chybu rekonstrukce.

Zdrojový kód

% Rekonstrukce

% Vyppracovali: T.Machacek M.Koblizek

im1=imread('k1.jpg');

im2=imread('k2.jpg');

u = [213.7757 248.2053 1

388.8416 183.7096 1

389.8417 114.4060 1

469.6757 84.6479 1

613.2651 348.1799 1

527.7550 384.6424 1

524.7146 440.0959 1

261.8330 368.9630 1]

uc = [135.7548 327.4878 1

245.4035 194.7429 1

241.2867 124.8000 1

288.9964 68.2094 1

540.4365 240.0000 1

499.6123 303.9771 1

495.2898 358.2171 1

228.0285 418.5920 1]

count = 8;

% ----------- zobrazeni bodu -----------

ax = [-5000, 1000, -2000, 2000];

close all;

figure(1);

imshow(im1);

hold on;

plot(u(:,1),u(:,2),'b*');

axis(ax);

figure(2);

imshow(im2);

hold on;

plot(uc(:,1),uc(:,2),'b*');

axis(ax);

% ----------- Vypocet epipolarni geometrie -----------

% vypocet matice A

A=[u(:,1).*uc(:,1), u(:,1).*uc(:,2), u(:,1), u(:,2).*uc(:,1), u(:,2).*uc(:,2), u(:,2), uc(:,1), uc(:,2), ones(size(u,1),1)];

% vypocet fundamentalni matice Q

[U,S,V] = svd(A);

q = null(U*S*V');

Q = reshape(q, 3, 3)';

% vypocet epipolar z fundamentalni matice matice

lambdax1=(Q*(uc'))';

lambdax2=u*Q;

% zobrazeni epipolar

figure(1);

hold on;

x1 = -5000;

x2 = 1000;

for i = 1:count

xx = [x1, x2];

yy = [-(lambda1(i,3)+lambda1(i,1)*x1)/lambda1(i,2), -(lambda1(i,3)+lambda1(i,1)*x2)/lambda1(i,2)];

plot(xx, yy); %zobrazeni usecky mezi body bod1, bod2

end;

figure(2);

hold on;

for i = 1:count

bod1 = [x1, x2];

bod2 = [-(lambda2(i,3)+lambda2(i,1)*x1)/lambda2(i,2), -(lambda2(i,3)+lambda2(i,1)*x2)/lambda2(i,2)];

plot(bod1, bod2);

end;

% vypocet epipolu

[U,S,V] = svd(Q);

e1 = U(:,3); % e1 jepoloha 2. kamery v 1.

e2 = V(:,3); % e2 jepoloha 1. kamery v 2.

%treti souradnice musi byt 1

e1 = e1/e1(3);

e2 = e2/e2(3);

%zobrazeni epipolu

figure(1);

hold on;

plot(e1(1), e1(2), 'r*');

figure(2);

hold on;

plot(e2(1), e2(2), 'r*');

% ----------------- Rekonstrukce z kalibrovanych kamer -----------------

% nacteni zadanych P1, P2

load pointcal_mat;

%priprava dat pro vypocet

A1 = P1(1:3, 1:3);

b1 = P1(:, 4);

A2 = P2(1:3, 1:3);

b2 = P2(:, 4);

A11 = inv(A1);

A22 = inv(A2);

Ab1 = A11 * b1;

Ab2 = A22 * b2;

X = [];

for i = 1:count

Au1 = A11 * u(i,:)';

Au2 = A22 * uc(i,:)';

pomL = [ Au1'*Au1, -Au2'*Au1; Au1'*Au2, -Au2'*Au2 ];

pomR = [ (Ab1'-Ab2')*Au1; (Ab1'-Ab2')*Au2 ];

aa = inv(pomL) * pomR;

X = [ X; (([aa(1) * Au1 - Ab1; 1]+[aa(2) * Au2 - Ab2; 1])/2)' ];

%presny bod je ziskan z prumeru rekonstrukce bodu z obou prumetu

end;

%vykresleni 3D rekonstruovanych bodu

figure(3);

plot3(X(:,1), X(:,2), X(:,3), 'b*');

grid on;

axis equal;

%reprojekce z rekonstruovanych bodu zpet do obou obrazku

newu = ( (P1*X') ./ ([1;1;1]*[P1(3,:) * X']) )';

newuc = ( (P2*X') ./ ([1;1;1]*[P2(3,:) * X']) )';

%zobrazeni reprojekce bodu do obrazku

figure(1);

hold on;

plot(newu(:,1),newu(:,2),'y*');

figure(2);

hold on;

plot(newuc(:,1),newuc(:,2),'y*');

% -------------- Projektivni rekonstrukce --------------

%volba P1 = [ I3 O3 ]

myP1 = [ 1 0 0 0;

0 1 0 0;

0 0 1 0 ];

%sestaveni matice W

f1 = Q(1, :);

f2 = Q(2, :);

f3 = Q(3, :);

f0 = zeros(1, 3); %nuly

W = [ f1 f0 f0;

f0 f2 f0;

f0 f0 f3;

f2 f1 f0;

f0 f3 f2;

f3 f0 f1 ];

%vypocet matice A

[ U, D, V ] = svd(W);

AA = V(:, 6:9);

%baze W (9x4)

%matici A volime jako linearnikombinaci vektoru baze W takovou, ktera je nejmene singularni

%budeme to zkouse 100x

mx=0;

for i=1:100

k = rand(4,1);

tmpA = AA * koef;

tmpA = reshape(tmpA, 3, 3);

D = svd(tmpA);

res = D(end)/D(1);

if ((res > mx) | (i==1))

mx = res;

myA = tmpA;

end;

%vypocet P2 = [ A b ] b=e2'

myP2 = [ myA e2' ];

%rekonstrukce bodu 3D sceny z bodu v obrazku pomoci nasich matic myP1, myP2

A1 = myP1(1:3, 1:3);

b1 = myP1(:, 4);

A2 = myP2(1:3, 1:3);

b2 = myP2(:, 4);

A11 = inv(A1);

A22 = inv(A2);

Ab1 = A11 * b1;

Ab2 = A22 * b2;

myX = []; %rekonstruovane body

for i = 1:count

Au1 = A11 * u(i,:)';

Au2 = A22 * uc(i,:)';

pomL = [ Au1'*Au1, -Au2'*Au1; Au1'*Au2, -Au2'*Au2 ];

pomR = [ (Ab1'-Ab2')*Au1; (Ab1'-Ab2')*Au2 ];

aa = inv(pomL) * pomR;

myX = [ myX; (([aa(1) * Au1 - Ab1; 1]+[aa(2) * Au2 - Ab2; 1])/2)' ];

end;

% reprojekce rekonstruovanych bodu zpet do obrazku

newu = ( (myP1*myX') ./ ([1;1;1]*[myP1(3,:) * myX']) )';

newuc = ( (myP2*myX') ./ ([1;1;1]*[myP2(3,:) * myX']) )';

%zobrazeni reprojekce bodu

figure(1);

hold on;

plot(newu(:,1),newu(:,2),'c*');

figure(2);

hold on;

plot(newuc(:,1),newuc(:,2),'c*');

% vypocet matice H

% sestaveni matice Z pro pet bodu podle prednasek

Z = [];

for i = 1:5

recX = X(i, :);

smyX = -myX(i, :)';

rZ = zeros(1,4);

sZ = zeros(4,1);

poml = [];

pomr = [];

%pro pet sloupcu leve

for j = 1:5

if (i == j)

poml = [poml, smyX];

else

poml = [poml, sZ];

end;

%pro ctyri ctyr-sloupce prave

for j = 1:4

sloupec = [];

%pro ctyri radky prave

for k = 1:4

if (j == k)

sloupec = [ sloupec; recX ];

else

sloupec = [ sloupec; rZ ];

end;

pomr = [pomr, sloupec];

end;

Z = [ Z; [poml, pomr] ];

end;

% vypocet mi-cek a matice H

h = null(Z); %odhad matice h

mis = h(1:5);

H = reshape(h(6:21), 4, 4)';

% prepocet nasi rekonstrukce do sceny pomoci matice H

recX = [];

for i = 1:size(mis)

newX = inv(H) * mis(i) * (myX(i,:)');

recX = [ recX; newX' ];

end;

%zobrazeni nasich rekonstruovanych bodu z matic myP1 a myP2 a H

figure(3);

hold on;

plot3(recX(:,1), recX(:,2), recX(:,3), 'g*');

grid on;

axis equal;

%reprojekce nasich rekonstruovanych bodu do obrazku

newu = ( (P1*recX') ./ ([1;1;1]*[P1(3,:) * recX']) )';

newuc = ( (P2*recX') ./ ([1;1;1]*[P2(3,:) * recX']) )';

%zobrazeni reprojekce bodu do obrazku

figure(1);

hold on;

plot(newu(:,1),newu(:,2),'m*'); % zobrazeni bodu

figure(2);

hold on;

plot(newuc(:,1),newuc(:,2),'m*'); % zobrazeni bodu