本練習程式是受到了這個老外博文的啟發，感覺挺有意思，就嘗試了一下。他用的是opencv，我這裡用的是matlab。

過去寫過透視投影，當時是用來做傾斜校正的，這次同樣用到了透視投影，不過更有意思，是將一張影象貼到另一張影象上。

兩個透視投影都需要先計算投影矩陣，傾斜校正那一篇是通過解線性方程組求的變換矩陣，而這一篇是通過奇異值分解求的變換矩陣。

為了對齊兩張影象，還需要對投影后的影象做一次仿射變換，其實就是座標平移。

這裡做投影和仿射直接呼叫了matlab的系統函式，方便一些。

還是先介紹下如何投影吧：

比如我們解單應性方程組或奇異值分解得到了投影矩陣：

那麼變換方程就可以寫為：

其中，x,y為原影象畫素座標，X,Y為目標影象畫素座標。

做一些正反點匹配變換的時候我會用這個公式，不過大多數情況下在做影象投影的時候，我並不用上面的公式，而是用下面這個公式：

其中，x,y為原影象畫素座標，X,Y為目標影象畫素座標。

第一個公式計算後圖像會出現空洞，第二個公式計算量明顯變大了，各有利弊吧。

這麼複雜的一個公式，我可推不出來，是用Mathematica推的，這個軟體真是極大的提高了我的效率。

話說matlab應該也能推，不過matlab的符號計算畢竟不是強項，我也就沒去了解matlab的這項功能。

看看效果吧。

原圖：

廣告牌：

結果：

matlab程式碼如下：

main.m

 1 clear all;close all;clc;
 
 2 
 3 img1=imread('lena.jpg');
 4 [h1 w1]=size(img1);
 5 mask=uint8(ones(h1,w1));    %二值模板，方便最後的合成
 6 
 7 img2=imread('pai.jpg');
 8 [h2 w2]=size(img2);
 9 
10 imshow(img1);
11 figure;imshow(img2);
12 
13 p1=[1,1;w1,1;1,h1;w1,h1];
14 p2=ginput();        %依次點選公告牌左上、右上、左下、右下
15 
16 T=calc_homography(p1,p2);   %計算單應性矩陣
 
17 T=maketform('projective',T);   %投影矩陣
18 
19 [imgn X Y]=imtransform(img1,T);     %投影
20 mask=imtransform(mask,T);
21 
22 T2=eye(3);
23 if X(1)>0, T2(3,1)= X(1); end
24 if Y(1)>0, T2(3,2)= Y(1); end
25 T2=maketform('affine',T2);      %仿射矩陣
26 
27 imgn=imtransform(imgn,T2,'XData',[1 w2],'YData',[1 h2]);    %仿射
28 mask=imtransform(mask,T2,'XData',[1 w2],'YData',[1 h2]);
29 
30 img=img2.*(1-mask)+imgn.*mask;  %合成
31 figure;imshow(img,[])

calc_homography

 1 function T = calc_homography(points1, points2)
 2 
 3     xaxb = points2(:,1) .* points1(:,1);
 4     xayb = points2(:,1) .* points1(:,2);
 5     yaxb = points2(:,2) .* points1(:,1);
 6     yayb = points2(:,2) .* points1(:,2);
 7 
 8     A = zeros(size(points1, 1)*2, 9);
 9     A(1:2:end,3) = 1;
10     A(2:2:end,6) = 1;
11     A(1:2:end,1:2) = points1;
12     A(2:2:end,4:5) = points1;
13     A(1:2:end,7) = -xaxb;
14     A(1:2:end,8) = -xayb;
15     A(2:2:end,7) = -yaxb;
16     A(2:2:end,8) = -yayb;
17     A(1:2:end,9) = -points2(:,1);
18     A(2:2:end,9) = -points2(:,2);
19 
20     [junk1,junk2,V] = svd(A);
21     h = V(:,9) ./ V(9,9);
22     T= reshape(h,3,3);
23 end