過去也寫過透視變換,當時算法真是弱爆了,我竟然會通過兩次變換。不過那引用的三篇文章都是非常好的文章,直到今天我才看明白。所謂的傾斜校正,一定要有標定點,將一個傾斜的矩形變為不傾斜的。因此可以從原四邊形四個點和新矩形四個點得到一個變換矩陣,根據這個矩陣再作用到全局圖像就可以了。詳細原理在這里,MIT的,我也不會比他介紹的更好了,還是看原版的好。
我這里的代碼完全就是按照MIT那篇文章的原理實現的,不過因為Matlab細節的原因,我把公式中x和y位置互換了:
clear all; close all; clc; img= imread('rect.bmp'); img= rgb2gray(img); imshow(mat2gray(img)); [M N] = size(img); dot=ginput(); %取四個點,依次是左上,右上,左下,右下,這里我取的是書的四個角 w=round(sqrt((dot(1,1)-dot(2,1))^2+(dot(1,2)-dot(2,2))^2)); %從原四邊形獲得新矩形寬 h=round(sqrt((dot(1,1)-dot(3,1))^2+(dot(1,2)-dot(3,2))^2)); %從原四邊形獲得新矩形高 y=[dot(1,1) dot(2,1) dot(3,1) dot(4,1)]; %四個原頂點 x=[dot(1,2) dot(2,2) dot(3,2) dot(4,2)]; %這里是新的頂點,我取的矩形,也可以做成其他的形狀 %大可以原圖像是矩形,新圖像是從dot中取得的點組成的任意四邊形.:) Y=[dot(1,1) dot(1,1) dot(1,1)+h dot(1,1)+h]; X=[dot(1,2) dot(1,2)+w dot(1,2) dot(1,2)+w]; B=[X(1) Y(1) X(2) Y(2) X(3) Y(3) X(4) Y(4)]'; %變換后的四個頂點,方程右邊的值 %聯立解方程組,方程的系數 A=[x(1) y(1) 1 0 0 0 -X(1)*x(1) -X(1)*y(1); 0 0 0 x(1) y(1) 1 -Y(1)*x(1) -Y(1)*y(1); x(2) y(2) 1 0 0 0 -X(2)*x(2) -X(2)*y(2); 0 0 0 x(2) y(2) 1 -Y(2)*x(2) -Y(2)*y(2); x(3) y(3) 1 0 0 0 -X(3)*x(3) -X(3)*y(3); 0 0 0 x(3) y(3) 1 -Y(3)*x(3) -Y(3)*y(3); x(4) y(4) 1 0 0 0 -X(4)*x(4) -X(4)*y(4); 0 0 0 x(4) y(4) 1 -Y(4)*x(4) -Y(4)*y(4)]; fa=inv(A)*B; %用四點求得的方程的解,也是全局變換系數 a=fa(1);b=fa(2);c=fa(3); d=fa(4);e=fa(5);f=fa(6); g=fa(7);h=fa(8); rot=[d e f; a b c; g h 1]; %公式中第一個數是x,Matlab第一個表示y,所以我矩陣1,2行互換了 pix1=rot*[1 1 1]'/(g*1+h*1+1); %變換后圖像左上點 pix2=rot*[1 N 1]'/(g*1+h*N+1); %變換后圖像右上點 pix3=rot*[M 1 1]'/(g*M+h*1+1); %變換后圖像左下點 pix4=rot*[M N 1]'/(g*M+h*N+1); %變換后圖像右下點 height=round(max([pix1(1) pix2(1) pix3(1) pix4(1)])-min([pix1(1) pix2(1) pix3(1) pix4(1)])); %變換后圖像的高度 width=round(max([pix1(2) pix2(2) pix3(2) pix4(2)])-min([pix1(2) pix2(2) pix3(2) pix4(2)])); %變換后圖像的寬度 imgn=zeros(height,width); delta_y=round(abs(min([pix1(1) pix2(1) pix3(1) pix4(1)]))); %取得y方向的負軸超出的偏移量 delta_x=round(abs(min([pix1(2) pix2(2) pix3(2) pix4(2)]))); %取得x方向的負軸超出的偏移量 inv_rot=inv(rot); for i = 1-delta_y:height-delta_y %從變換圖像中反向尋找原圖像的點,以免出現空洞,和旋轉放大原理一樣 for j = 1-delta_x:width-delta_x pix=inv_rot*[i j 1]'; %求原圖像中坐標,因為[YW XW W]=fa*[y x 1],所以這里求的是[YW XW W],W=gy+hx+1; pix=inv([g*pix(1)-1 h*pix(1);g*pix(2) h*pix(2)-1])*[-pix(1) -pix(2)]'; %相當於解[pix(1)*(gy+hx+1) pix(2)*(gy+hx+1)]=[y x],這樣一個方程,求y和x,最后pix=[y x]; if pix(1)>=0.5 && pix(2)>=0.5 && pix(1)<=M && pix(2)<=N imgn(i+delta_y,j+delta_x)=img(round(pix(1)),round(pix(2))); %最鄰近插值,也可以用雙線性或雙立方插值 end end end figure; imshow(uint8(imgn));
程序效果:
原圖,這是本不錯的書
傾斜校正后
將來說不定結合sift算子和霍夫變換就能自動校正呢。
注:博客園有Bug,當你點代碼左下角復制代碼的時候,矩陣A的2,4,6,8行第一列會少復制一個0.
另一篇關於投影的文章見這里。