python计算机视觉课程实验_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-8

import numpy as np

from scipy import linalg

class 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.c

def 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 R

def 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).T

def 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 cam

import numpy as np

import pylab as pl

from PIL import Image

im = 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.T

P = 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-8

import homography

import camera

import sift

''' generate sift features'''

'''sift tmp.pgm --output=im1.sift --edge-thresh 10 --peak-thresh 5'''

sift.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=0

box_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])]).T

cam2.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]).T

t = 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()