python計算機視覺第四章

4.1針孔照相機模型

針孔照相機模型是計算機視覺中廣泛使用地照相機模型

其中P為照相機矩陣，維度為3×4

，其中R是描述照相機方向的旋轉矩陣，t是描述照相機中心位置的三維平移向量

K = np.array([[af,s,cx],[0,f,cy],[0,0,1]]) f為圖像平面和照相機中心間的距離（焦距）大多數情況s=0，a=1，cx與cy分別為圖像高度和寬度的一半

python代碼以及解釋

創建針孔照相機類

圖4-2代碼

def

#coding=utf-8import numpy as npfrom scipy import linalgclass camera(object): def __init__(self,P): self.P = P #照相機矩陣 self.K = None #內標定矩陣 self.R = None #旋轉矩陣 self.t = None #平移矩陣 self.c = None #照相機中心 def project(self,X): X進行投影，並進行坐標歸一化 x = np.dot(self.P,X) for i in range(3): x[i] /=x[2] return x def factor(self): #分解前3×3部分 K,R = linalg.rq(self.P[:,:3]) T = np.diag(np.sign(np.diag(K))) if np.linalg.det(T)<0: T[1,1]*=-1 self.K = np.dot(K,T) self.R = np.dot(T,R) self.t = np.dot(linalg.inv(self.K),self.P[:,3]) return self.K,self.R,self.t def center(self): 計算並返回照相機的中心 if self.c is not None: return self.c else: self.factor() self.c = -np.dot(self.R.T,self.t) return self.cdef rotation_matrix(a): R = np.eye(4) R[:3,:3] = linalg.expm([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]]) return Rdef cube_points(c,wid): p=[] p.append([c[0]-wid,c[1]-wid,c[2]-wid]) p.append([c[0]-wid,c[1]+wid,c[2]-wid]) p.append([c[0]+wid,c[1]+wid,c[2]-wid]) p.append([c[0]+wid,c[1]-wid,c[2]-wid]) p.append([c[0]-wid,c[1]-wid,c[2]-wid]) p.append([c[0]-wid,c[1]-wid,c[2]+wid]) p.append([c[0]-wid,c[1]+wid,c[2]+wid]) p.append([c[0]+wid,c[1]+wid,c[2]+wid]) p.append([c[0]+wid,c[1]-wid,c[2]+wid]) p.append([c[0]-wid,c[1]-wid,c[2]+wid]) p.append([c[0]-wid,c[1]-wid,c[2]+wid]) p.append([c[0]-wid,c[1]+wid,c[2]+wid]) p.append([c[0]-wid,c[1]+wid,c[2]-wid]) p.append([c[0]+wid,c[1]+wid,c[2]-wid]) p.append([c[0]+wid,c[1]+wid,c[2]+wid]) p.append([c[0]+wid,c[1]-wid,c[2]+wid]) p.append([c[0]+wid,c[1]-wid,c[2]-wid]) return np.array(p).Tdef my_calibration(sz): row,col = sz fx = 2555*col/2592 fy = 2586*row/1936 K = np.diag([fx,ft,1]) K[0,2] = 0.5*col K[1,2] = 0.5*row return K

import camera as camimport numpy as npimport pylab as plfrom PIL import Imageim = np.array(Image.open(../data/Model House/house.001.pgm).convert(L))points = np.loadtxt(../data/Model House/3D/house.p3d)points = np.hstack((points,np.ones((points.shape[0],1))))points = points.TP = np.hstack((np.eye(3),np.array(([0],[0],[-10]))))r = 0.05*np.random.rand(3)rot = cam.rotation_matrix(r)cam = cam.camera(P)x = cam.project(points)pl.figure()pl.gray()pl.subplot(1,3,1)pl.imshow(im)pl.subplot(1,3,2)pl.plot(x[0],x[1],*)pl.subplot(1,3,3)for t in range(20): cam.P = np.dot(cam.P,rot) x = cam.project(points) pl.plot(x[0],x[1],+) # 注意對於旋轉，t向量是保持不變的pl.show()

4.1.3 照相機矩陣的分解

對於P，恢復內參數K以及照相機的位置t和姿勢R。矩陣分塊操作稱為因子分解。使用RQ因子分解方法

且

根據要求，需要將K中對角元素設置為大於0，因此在使用linalg.rq函數以後，再使用diag函數獲得K對角元素構成的向量，隨後通過sign函數返回三個對角元素符號判斷（大於0為1，否則為0），再調用diag生成一個對角矩陣；最終使用det計算行列式的值從而完成判斷

4.1.4計算照相機中心

照相機中心c是一個三維點，滿足約束

#coding=UTF-8import homographyimport cameraimport sift generate sift featuressift tmp.pgm --output=im1.sift --edge-thresh 10 --peak-thresh 5sift.process_image(book_frontal.JPG,im0.sift) #生成sift特徵文本#lo 生成N*4矩陣 （x坐標，y坐標，尺度，主方向）lo,do = sift.read_features_from_file(im0.sift)#從文件中讀取特徵點坐標以及描述符向量sift.process_image(book_perspective.JPG,im1.sift) #生成sift特徵文本l1,d1 = sift.read_features_from_file(im1.sift)#從文本中讀取特徵點坐標以及描述符向量matches = sift.match_twosided(d0,d1)ndx = matches.nonzero()[0]fp = homnography.make_homog(lo[ndx,:2].T) #生成齊次坐標ndx2 = [int(matches[i) for i in ndx]tp = homnography.make_homog(l1[ndx,:2].T) #生成齊次坐標 就是在底層增加1×N的ones向量model = homography.RansacModel()H = hmnnography.H_from_ransc(fp,tp,model)[0] #得到fp to tp的單應性矩陣K = my_calibration((747,1000))box = cube_points([0,0,0.1],0.1)cam1 = camera.Camera(np.hstack(K,np.dot(K,np.array([[0],[0],[-1]]))))#type:camera#將底部點進行投影 將5個點首先齊次化，然後調用project進行點× 且z=0box_cam1 = camera.project(homography.make_homog(box[:,:5]))#使用H，將box_cam1的點變換到第二幅圖像中 使用normalize歸一化，變為齊次形式box_trans = homography.normalize(np.dot(H,box_cam1))#根據單應性矩陣，得到新cam的照相機矩陣 P=HP#由於前三列乘以K的逆，得到的是旋轉矩陣R，對於一般的旋轉矩陣，第三列與其他兩列呈現的是正交關係，因此將第三列替換為第一第二列的叉乘cam2 = camera.Camera(np.dot(H,cam1.P))A = np.dot(linalg.inv(K),cam2.P[:,:3])A = np.array([A[:,0],A[:,1],np.cross(A[:,0],A[:,1])]).Tcam2.P[:,:3] = np.dot(K,A)#將box進行投影（使用第二個照相機矩陣）box_cam2 = cam2.project(homography.make_homnog(box))#對於在z面上一個點 point 首先計算通過單應性矩陣得到的 然後使用cam2計算point 兩者進行比較point = np.array([1,1,0,1]).Tt = np.dot(H,cam1.project(point))print t/t[-1] #齊次坐標歸一化print cam2.project(point)

輸出為：

[ 1.47429285e+03 1.01011564e+02 1.00000000e+00]

[ 1.47429285e+03 1.01011564e+02 1.00000000e+00]

圖4-4

im0 = np.array(Image.open(book_frontal.JPG))im1 = np.array(Image.open(book_perspective.JPG))box = camera.cube_points([-0.5,-0.4,0.05],0.05)pl.figure()pl.gray()pl.subplot(2,2,1)pl.imshow(im0)pl.plot(box_cam1[0,:],box_cam1[1,:],linewidth=3)pl.subplot(2,2,2)pl.imshow(im1)pl.plot(box_trans[0,:],box_trans[1,:],linewidth=3)pl.subplot(2,1,2)pl.imshow(im1)pl.plot(box_cam2[0,:],box_cam2[1,:],linewidth=3)pl.show()