Počítačové vidění pro informatiku - Úloha 1

Josef Šivic

Contents

• Contents
• Zadání
• Odvození
• Řešení v matlabu
• Výsledky
• Bibliography

Zadání

Sestavte do mozaiky dva obrazy téže scény, jsou-li sejmuté ze stejného místa. V obrazech máte dány 4 korespondující body.

   
Odvození

Odvození vychází z [1] str. 450-453. Vyjdeme z rovnice (9.6) na str. 452.

$$z_c \left(\begin{array}{c} u\\ v\\ 1 \end{array}\right) =KR(X_w-C_w)$$ (1)

kde :
	K	je	kalibrační matice kamery
	R	je	rotační matice vyjadřující přechod od světových souřadnic
			k souřadnicím kamery
	Xw	je	bod vyjádřený ve světových souřadnicích
	Cw	je	počátek souřadného systému kamery vyjádřený ve
			světových souřadnicích
	(u,v,1)T	je	bod $\bf {\tilde{u}}$ ležící v obrazové rovině, vyjádřený v homogenních souřadnicích

Přepíšeme-li bod X_w do homogenních souřadnic, lze pak posunutí mezi počátkem světového systému souřadnic a systému souřadnic kamery vyjádřit pomocí maticového násobení. Rovnici (1) lze pak zapsat jako:

$$\alpha \left(\begin{array}{c} u\\ v\\ 1 \end{array}\right)... ...es1}) \left(\begin{array}{c} \bf {X_w}\\ 1 \end{array}\right)$$ (2)

kde :
	$(u,v,1)^T=\bf {\tilde{u}}$	je bod ležící v obrazové rovině, vyjádřený
		v homogenních souřadnicích
	$({\bf X_w}, 1)^T=(x,y,z,1)^T=\bf {\tilde{X}_w}$	je bod světových souřadnic, vyjádřený
		v homogenních souřadnicích

Máme-li tedy dva obrázky scény sejmuté různými kamerami z různých míst a zadáno Nkorespondujících bodů u_i^j ,kde j=1,2, $i=1 \ldots N$ platí rovnice

 

$$\alpha_i^1 \, \tilde{u}_i^1 K^1R^1(I\vert-{\bf C_w}^1) \tilde{X}_i$$ = (3)

$$\alpha_i^2 \, \tilde{u}_i^2 K^2R^2(I\vert-{\bf C_w}^2) \tilde{X}_i.$$ =

Sejmeme-li podle zadání úlohy obrázky ze stejného místa (ale s možným různým natočením a nakloněním) a ztotožníme-li počátek světového systému souřadnic a systému souřadnic kamery, rovnice

(3) se změní takto

 

$$\alpha_i^1 \, \tilde{u}_i^1 K^1R^1(I\vert{\bf0}) \tilde{X}_i$$ = (4)

$$\alpha_i^2 \, \tilde{u}_i^2 K^2R^2(I\vert{\bf0}) \tilde{X}_i.$$ =

Rovnice (4) lze pak dále upravovat

  $$\begin{eqnarray*} \alpha_i^1 \, \tilde{u}_i^1 &= & K^1R^1(I\vert{\bf0}) \left(... ...rray}{c} x_i/z_i\\ y_i/z_i\\ z_i/z_i\\ \end{array}\right)z_i,\\ \end{eqnarray*}$$

což lze přepsat jako

 

$$\alpha_i^1 \, \tilde{u}_i^1 K^1R^1 \tilde{x}_i \, z_i$$ = (5)

$$\alpha_i^2 \, \tilde{u}_i^2 K^2R^2 \tilde{x}_i \, z_i.$$ =

Vyjádříme-li z druhé rovnice z (5) $\tilde{x}_i.z_i$, dosadíme výsledek do první rovnice z (5) dostaneme

$$\frac{\alpha_i^1}{\alpha_i^2} \, \tilde{u}_i^1 = K^1R^1(R^2)^{-1}(K^2)^{-1} \, \tilde{u}_i^2. \\$$ (6)

Po zavedení $\gamma_i=\alpha_i^1/ \alpha_i^2$ a H=K¹R¹(R²)^-1(K²)^-1dostaneme hledaný vztah mezi korespondujícími body ve dvou obrazech sejmutých z jednoho místa

$$\gamma_i \, \tilde{u}_i^1 = H \, \tilde{u}_i^2. \\$$ (7)

Úkolem je tedy nalézt zobrazení $H=\{h_{ij}\}$.

Rovnici (7) lze přepsat jako

$$\alpha_i \left(\begin{array}{c} u_i\\ v_i\\ 1\end{array}\right)=H \left(\begin{array}{c} x_i\\ y_i\\ 1\end{array}\right).$$ (8)

Po úpravách rovnice (8 dostaneme pro každý bod $i=1 \ldots N$ dvojici rovnic

$${\tiny \begin{array}{cccccccccl} -x_i h_{11} & -y_i h_{12} ... ...v_i h_{31} & +y_i v_i h_{32} & +v_i h_{33} & = 0. \end{array}}$$ (9)

Takovouto soustavu lze převést do maticového tvaru

$$\underbrace{A}_{2N\times9} \cdot \underbrace{\left(\begin{ar... ...h_{12}\\ \vdots\\ h_{33}\end{array}\right)}_{9\times1}={\bf0}.$$ (10)

Aby byla soustava (10) řešitelná, potřebujeme znát alespoň 4 korespondující body. Známe-li právě 4 korespondující body, hod A=8 a soustava má právě jedno netriviální jednodimenzionální řešení $\lambda {\bf h}, \lambda \not= 0$.

Máme-li více korespondujících bodů, hod A může být 9 a v takovém případě má soustava (10) pouze triviální řešení. Cílem tedy je vhodným způsobem snížit hod A na 8. Jedna z možností je rozložit A pomocí SVD na A=UDV^T. Vynulováním nejmenšího singulárního čísla dostaneme matici

$$\bar{A}=U\left(\begin{array}{cccc} \sigma_1\\ &\ddots&\\ & & \sigma_{n-1}\\ & & & 0\end{array}\right)V^T.$$

kde $hod\bar{A}=8$ a $\bar{A}$ je k A 'nejbližší' ve smyslu frobeniovy normy viz [1] na str. 456.

Po dosazení $\bar{A}$ místo A do (10) dostaneme

$$U\left(\begin{array}{cccc} \sigma_1\\ &\ddots&\\ & & \sigma_... ... v}_8^T\\ {\bf v}_9^T\end{array}\right) \cdot {\bf h} = {\bf0}$$ (11)

Hledáme tedy takový vektor ${\bf h}$, který je ortogonální k ${\bf v}_1 \ldots {\bf v}_8$, tj. musí platit ${\bf v}_i^T{\bf h}=0$ pro $i=1 \ldots 8$. Na ortogonalitě k ${\bf v}_9$nezáleží, protože tam je nulovost rovnice zajištěna ( $\sigma_9=0$). Je vidět, že řešením je vektor ${\bf v}_9$.

Transformační matice má tedy tvar

$$H=\left(\begin{array}{ccc} v_{9,1} & v_{9,2} & v_{9,3} \\ v... ...{9,6} \\ v_{9,7} & v_{9,8} & v_{9,9} \\ \end{array}\right).$$

Řešení v matlabu

% PVI Uloha 1 - Mozaika
%
% Vstup:  2 obrazky sceny sejmute z jednoho mista:     im1, im2
%         N korespondujicich bodu:                     u1 , u2
%
% Vystup: Vysledny obrazek obsahujici slozene obrazky: im
%
% Autor:  Josef Sivic
% Datum:  4/4/2000 


%nacti obrazky
im1=double(rgb2gray(imread('panorama1.jpg')));
im2=double(rgb2gray(imread('panorama2.jpg')));

im1=im1/max(im1(:));
im2=im2/max(im2(:));

%nacti korespondujici body
load panoramap.mat

%ukaz obrazky a body
imshow(im1)
hold on
plot(I1(:,1),I2(:,1), '*');

figure
imshow(im2)
hold on
plot(I1(:,2),I2(:,2), '*');


%vytvor velky obrazek a uloz do nej prvni obrazek
im=zeros(size(im1)*2);
im( size(im,1)/4+(1:size(im1,1)) ,(size(im,2)-size(im1,2))+1:size(im,2))=im1(:,:);
figure
imshow(im);


%transoformuj souradnice u, do obrazku im
u1=[I2(:,1)' ; I1(:,1)']; %im1
u2=[I2(:,2)' ; I1(:,2)']; %im2 

o2=[size(im,1)/4+1,size(im,2)-size(im1,2)+1]; % pocatek im1 v im
u(1,:)=(u1(1,:)+o2(1));
u(2,:)=(u1(2,:)+o2(2));

hold on
plot(u(2,:),u(1,:),'b*')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% matice homografie H: im2 -> im     %
%                      u2  -> u      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if 0
 %vytvori matici A 
 N=size(u2,2);
 A=zeros(8,2*N);
 for i=1:N
    A(2*i-1,:)  =[u1(1,i) u1(2,i) 1 0         0       0  
                 -u1(1,i)*u2(1,i)  -u1(2,i)*u2(1,i)  -u2(1,i)];
    A(2*i,:)=    [0       0       0     u1(1,i) u1(2,i) 1  
                  -u1(1,i)*u2(2,i)  -u1(2,i)*u2(2,i)  -u2(2,i)]; 
 end;
 
 [V,D]=eig(A'*A);
 
end


H=uu2homog(u2,u,'norm');
H_=inv(H);



%Transformuj rohy im2 do im v poradi UL, LL, LR, UP
%1.radek radkova souradnic
%2.radek sloupcova souradnice
im2_corn=[1  1; size(im2,1) 1; size(im2,1) size(im2,2);  1 size(im2,2)]';

%transformuj rohy im2 do souradnic im
im2_corn_t=p2e(H * e2p(im2_corn));
plot(im2_corn_t(2,:),im2_corn_t(1,:),'m*');


%opsany obdelnik uvnitr ktereho se pracuje
mir=floor(min(im2_corn_t(1,:)));
mar=ceil(max(im2_corn_t(1,:)));
mic=floor(min(im2_corn_t(2,:)));
mac=ceil(max(im2_corn_t(2,:)));

im2_rect=[mir mic;  mar mic;  mar mac;  mir mac]';
plot(im2_rect(2,[1:4 1]),im2_rect(1,[1:4 1]),':r');

rect=im2_rect;       



%%%%%%%%%%%%%%%%%%%%
% Maticovy vypocet %
%%%%%%%%%%%%%%%%%%%%

% vytvor index idx
% pocet sloupcu je stejny jako pocet prvku v opsanem obdelniku
% pocet radku je 2:
% 1. radek je radkovy index prvku opsaneho obdelniku 
% 2. radek je sloupcovy index  ...

% nezobrazuj posledni sloupce (v obrazku chybi)
mac=mac-10;

% napln pole idx
[rows,cols]=meshgrid(mir:mar, mic:mac);
idx=[rows(:)'; cols(:)'];

% transformace souradnic im do souradnic im2
idx_t=p2e(H_*e2p(idx));

% vybrani jen platnych souradnic, takovych, ktere opravdu lezi v im2
% opsany obdelik pokryva vice
valid=((idx_t(1,:)>=2) & (idx_t(2,:)>=2) & 
       (idx_t(1,:)<=size(im2,1)) & (idx_t(2,:)<=size(im2,2)));



%%%%%%%%%%%%%%%%%%%% 
% prepis im2 do im %
%%%%%%%%%%%%%%%%%%%%

% v idx a idx_t jsou ulozeny souradnice sobe odpovidajicich pixelu v im a im2
% Chceme-li rychle pixely z im2 do im pridat je potreba ze souradnic ve 
% forme A(c,r) udelat souradnice a(c*m+r)

% mam-li prvek v matici A na miste A(r,c) (matice ma rozmery (m,n))
% Necht a=A(:), pak A(r,c) je na miste a(c*m + r)

ix=idx(2,:)*size(im,1) + idx(1,:);
ix_t=round(idx_t(2,:)-1)*size(im2,1) + round(idx_t(1,:));

%prepis obrazku
im(ix(valid))=im2(ix_t(valid));




if 0
%vypocet for cyklem
for i=mir:mar      %pro vsechny radky z obdelniku v im
   for j=mic:mac   %pro vsechny sloupce z obdelniku v im
     
      point=p2e(H_*e2p([i j]'));
      r=round(point(1));
      c=round(point(2)); 
      if (r<=size(im2,1)) & (c<=size(im2,2) & r>=1 & c>=1)
         im(i,j)=im2(r,c);
      end   
   end  
end

end

% zobraz vysledek
figure
hold on
imshow(im);
axis ij;

Výsledky

Obrázky na sobě ve sýsledném obrázku 3 úplně přesně nesedí. Chyba vznikla nejspíše z nepřesně určených korespondujících bodů. Určíme-li na vstupních datech co nejpřesněji více korespondujících bodů, výsledky by měly být lepší. Problém algoritmu ale je, že i jen jeden hodně nepřesně zadaný bod negativně ovlivní celý výpočet matice H.

  

Figure 1: Vstupní obrázky s vyznačenými korespondencemi im1(vlevo), im2(vpravo)

  

Figure 2: Výstupní obrázek im obsahující vstupní obrázek im1, vyznacene transformovane rohy obrázku im2(*) a jim opsaný obdélník

  

Figure 3: Výsledný obraz obsahující vstupní obraz im1 (vpravo) a transformovaný obraz im2 (vlevo). Byly použity 4 korespondující body.

Bibliography

1

©onka M., Hlavaè V., Boyle R.,
Image Processing, Analysis, and Machine Vision,
PWS, Boston, USA, 1998