ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ

fakulta elektrotechnická

Semestrální práce z předmětu

Vypracovali: David Vilímek

Martin Čadík

Vždy napíši nejprve teorii poté odpovídající kód a poté obrázky.

u1, u2..jsou vektory bodů v prvním a druhém pohledu

S(T), R.. jsou další matice určující postavení kamer

u1^T*K1^-T*S(T)*R^-1*K2^-1*u2=0

u1^T*Q*u2=0

u1^T*l 1=0, l 2*u2=0

(u1_i, v1_i, 1)*Q*(u2_i, v2_i, 1)^T=0

(u1_i* u2_i, u1_i* v2_i, u1_i, v1_i* u2_i, v1_i* v2_i, v1_i, u2_i, v2_i, 1)* Q11

Q12

Stačí tedy 8 bodů.

Epipoláry

(u1_i, v1_i, 1)* l 1=0

u1_i* l 1¹+ v1_i *l 1²+ l 1³=0

x*a+y*b+c=0

load kostky2.cor -mat

a = CORR_EXCHANGE.d.corr;

u = [a(:,1,1) a(:,1,2) ones(12,1)]; %prvni obraz x,y - j, i

uc = [a(:,2,1) a(:,2,2) ones(12,1)]; %druhy obraz x,y - j, i

figure;

im1=(imread('k1.jpg'));

imshow(im1);

axis([-2000 2000 -2000 2000]);

hold;

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

hold off;

figure;

im2=(imread('k2.jpg'));

imshow(im2);

axis([-1000 1000 -1000 1000]);

hold;

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

hold off;

count = size(u,1);

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)];

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

S(9,9) = 0;

A2 = U*S*V';

% vypocet Q

q = null(A2);

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

[ux,s,v] = svd(Q);

s(3,3) = 0;

Qsvd = ux*s*v';

u1=[u(:,2) u(:,1) ones(size(u,1),1)];

uc1=[uc(:,2) uc(:,1) ones(size(uc,1),1)];

lambda1 = []; %v prvnim obrazku

lambda2 = []; %v druhem obrazku

for i = 1:count

lambda1 = [lambda1; (Q*uc(i,:)')'];

lambda2 = [lambda2; u(i,:)*Q];

end;

b1 = [];

b2 = [];

k = 5000;

for i = 1 : count

b1 = [b1 [k; -k]]; %x1 x2

b2 = [b2 [-(lambda1(i,3)+k*lambda1(i,1))/lambda1(i,2); -(lambda1(i,3)-k*lambda1(i,1))/lambda1(i,2)]]; %y1 y2

end

figure(1);

hold on;

axis([0 2600 0 700]);

axis on;

plot(b1, b2);

e1 = null(Qsvd');

e1 = e1/e1(3,1);

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

b1 = [];

b2 = [];

k = 5000;

for i = 1 : count

b1 = [b1 [k; -k]]; %x1 x2

b2 = [b2 [(-lambda2(i,3)-k*lambda2(i,1))/lambda2(i,2); (-lambda2(i,3)+k*lambda2(i,1))/lambda2(i,2)]]; %y1 y2

end

figure(2);

hold on;

axis([0 1600 -100 600]);

axis on;

plot(b1, b2);

e2 = null(Qsvd);

e2 = e2/e2(3,1);

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

Nyní reprojekce "do světových" souřadnic se zadaných matic pro kamery a zobrazení 3D

a 1*u1=P1*X

a 2*u2=P2*X

My počítáme X

load pointcal_mat;

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

b1 = P1(:, 4);

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

b2 = P2(:, 4);

A11 = inv(A1);

A22 = inv(A2);

u1=[u1(:,2) u1(:,1) u1(:,3)];

uc1=[uc1(:,2) uc1(:,1) uc1(:,3)];

X1 = [];

X2 = [];

for i = 1:count

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

Ab1 = A11 * b1;

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

Ab2 = A22 * b2;

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

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

alphas = inv(LEFT) * RIGHT;

X1 = [ X1; [alphas(1) * Au1 - Ab1; 1]' ];

X2 = [ X2; [alphas(2) * Au2 - Ab2; 1]' ];

end;

%presny bod jako prumer obou prumetu

X = (X1 + X2) / 2;

figure(3);

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

grid on;

axis equal;

%normalizace na 3.slozku == 1

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

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

%zobrazeni pretransformovanych bodu

%pause;

figure(1);

hold on;

plot(newPoints1(:,1),newPoints1(:,2),'y*'); % zobrazeni bodu

figure(2);

hold on;

plot(newPoints2(:,1),newPoints2(:,2),'y*'); % zobrazeni bodu

u1^T*Q*u2=0

(P1*X_i)^T*Q*(P2*X_i)=0

X_i ^T *P1 ^T *Q*P2*X_i=0

a 1_i*u1_i=P1*X_i=(I3|O3)*X_i

a 2_i*u2_i=P2*X_i

X_i ^T *(I3|O3)^T *Q*P2*X_i=0

X_i =(x ^T x ) ^T , P2=(A | b), b =e2

(x ^T x ) *(I3|O3)^T *Q*(A | b)* (x ^T x ) ^T =0

x ^T*Q*A* x + x ^T*Q*b*x =0

x ^T*Q*A* x =0

pak Q*A je antisymetrická

q1

q2

q3 a1 a1

q2 q1 a2 =0, W* a2 =0, hod(W)=5

q3 q2 a3 a3

q3 q1

%volba P1 = [ I3 O3 ]

myP1 = [1 0 0 0;

0 1 0 0;

0 0 1 0 ];

%vypocet P2 = [ A b ] ---

myb = e2; % e2 = P1 * c1 = b

%sestaveni matice W

f1 = Q(1, :);

f2 = Q(2, :);

f3 = Q(3, :);

f0 = zeros(1, 3); %same nuly

W = [ f1 f0 f0;

f0 f2 f0;

f0 f0 f3;

f2 f1 f0;

f0 f3 f2;

f3 f0 f1 ];

%vypocet pomocne matice A pomoci SVD z W

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

pomA = V(:, 6:9); %baze W (9x4), matice A (jako radkovy vektor) je LK radkovych vektoru teto baze

%vypocet nejvhodnejsich koeficientu (alfa) pro zisk A - chceme co nejmene singularni => male rozdily v singul. cislech

%zkousime maxRuns-krat nahodne volit koeficienty, bereme nejlepsi

maxRuns = 1000; %kolikrat zkousim

koef = rand(4,1); %nahodne koeficienty (4x1)

testA = pomA * koef; %dosud nejlepsi matice

myA = reshape(testA, 3, 3);

D = svd(myA); %!vraci vektor singul. cisel (diagonalu D)

bestRegul = D(end)/D(1); %dosud nejlepsi hodnota singularity (chceme co nejblizsi nejmensi a nejvetsi cislo-hodnoty na diag klesaji

for i=1:maxRuns

koef = rand(4,1);

testA = pomA * koef;

testA = reshape(testA, 3, 3);

D = svd(testA);

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

if (pomRegul > bestRegul)

bestRegul = pomRegul;

myA = testA;

end;

end;

%sestaveni P2

myP2 = [ myA myb ];

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

b1 = myP1(:, 4);

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

b2 = myP2(:, 4);

A11 = inv(A1);

A22 = inv(A2);

%vypocet bodu "ve svete" z rovnice - viz papir (bez nej ani ranu)

X1 = [];

X2 = [];

for i = 1:count

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

Ab1 = A11 * b1;

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

Ab2 = A22 * b2;

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

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

alphas = inv(LEFT) * RIGHT;

X1 = [ X1; [alphas(1) * Au1 - Ab1; 1]' ];

X2 = [ X2; [alphas(2) * Au2 - Ab2; 1]' ];

end;

%presny bod jako prumer obou prumetu

myX = (X1 + X2) / 2;

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

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

%zobrazeni pretransformovanych bodu

%pause;

figure(1);

hold on;

plot(newPoints1(:,1),newPoints1(:,2),'r*'); % zobrazeni bodu

figure(2);

hold on;

plot(newPoints2(:,1),newPoints2(:,2),'r*'); % zobrazeni bodu

m _i*X1_i=H*X2_i

X = X(2:6,:);

myX = myX(2:6,:);

Z = [];

move = 0;

%pro pet bodu

for i = 1:5

bod = move + i;

rX = X(bod, :); %bod ve scene - radkovy vektor

smyX = -myX(bod, :)';%bod v rekonstrukci - sloupcovy vektor

rZ = zeros(1,4); %nuly v radku

sZ = zeros(4,1); %nuly ve sloupci

left = [];

right = [];

%pro pet sloupcu leve

for j = 1:5

if (i == j)

left = [left, smyX];

else

left = [left, sZ];

end;

end;

%pro ctyri ctyr-sloupce prave

for j = 1:4

columnRight = [];

%pro ctyri radky prave

for k = 1:4

if (j == k)

columnRight = [ columnRight; rX ];

else

columnRight = [ columnRight; rZ ];

end;

end;

right = [right, columnRight];

end;

Z = [ Z; [left, right] ];

end;

% vypocet mi-cek a matice H

test = null(Z);

mis = test(1:5);

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

% prepocet nasi rekonstrukce do sceny pomoci H matice ---------------------------------

%vypocet realXi = inv(H) * MIi * myXi pro body, pro ktere mame MIi

myRealX = [];

for i = 1:size(mis)

bod = move + i;

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

myRealX = [ myRealX; newX' ];

end;

%zobrazeni vypoctenych bodu

figure(3);

hold on;

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

% zobrazeni nasich bodu v realnem svete (transformaci pres H) zpet realnymi kamerami ----------

% prepocet bodu (alfa*u = P*X1) s normalizaci na 3.slozku == 1

newPoints1 = ( (P1*myRealX') ./ ([1;1;1]*[P1(3,:) * myRealX']) )';

newPoints2 = ( (P2*myRealX') ./ ([1;1;1]*[P2(3,:) * myRealX']) )';

%zobrazeni pretransformovanych bodu

figure(1);

hold on;

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

figure(2);

hold on;

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

Proto jsme poté zkoumali závislost odchylek na volbě vstupních bodů. Je zřejmé, že odchylka velice roste v případě, že mezi vstupními body je takový, jehož souřadnice jsou zadány s hrubou nepřesností, a to i pro ostatní body zadané relativně přesně.