張正友相機標定法是張正友教授1998年提出的單平面棋盤格的相機標定方法。傳統(tǒng)標定法的標定板是需要三維的，需要非常精確，這很難制作，而張正友教授提出的方法介于傳統(tǒng)標定法和自標定法之間，但克服了傳統(tǒng)標定法需要的高精度標定物的缺點，而僅需使用一個打印出來的棋盤格就可以。同時也相對于自標定而言，提高了精度，便于操作。因此張氏標定法被廣泛應(yīng)用于計算機視覺方面。

傳統(tǒng)標定法的的標定板

張正友標定法的標定板

今天，我們就來講解一下張氏標定法的原理和實現(xiàn)，學(xué)會之后，我們就可以自己去制作一個棋盤標定板，然后拍照，標定自己手機相機的參數(shù)啦！

一、相機標定介紹

二、算法原理

1.整體流程

2.模型假設(shè)

3.模型求解

（1）內(nèi)外參數(shù)求解

（2）畸變系數(shù)求解

（3）精度優(yōu)化

三、算法實現(xiàn)

1.main.py

2.homography.py

4.extrinsics.py

5.distortion.py

6.refine_all.py

7.結(jié)果

一、相機標定介紹

相機標定指建立相機圖像像素位置與場景點位置之間的關(guān)系，根據(jù)相機成像模型，由特征點在圖像中坐標與世界坐標的對應(yīng)關(guān)系，求解相機模型的參數(shù)。相機需要標定的模型參數(shù)包括內(nèi)部參數(shù)和外部參數(shù)。

針孔相機成像原理其實就是利用投影將真實的三維世界坐標轉(zhuǎn)換到二維的相機坐標上去，其模型示意圖如下圖所示：

從圖中我們可以看出，在世界坐標中的一條直線上的點在相機上只呈現(xiàn)出了一個點，其中發(fā)生了非常大的變化，同時也損失和很多重要的信息，這正是我們3D重建、目標檢測與識別領(lǐng)域的重點和難點。實際中，鏡頭并非理想的透視成像，帶有不同程度的畸變。理論上鏡頭的畸變包括徑向畸變和切向畸變，切向畸變影響較小，通常只考慮徑向畸變。

徑向畸變：徑向畸變主要由鏡頭徑向曲率產(chǎn)生（光線在遠離透鏡中心的地方比靠近中心的地方更加彎曲）。導(dǎo)致真實成像點向內(nèi)或向外偏離理想成像點。其中畸變像點相對于理想像點沿徑向向外偏移，遠離中心的，稱為枕形畸變；徑向畸點相對于理想點沿徑向向中心靠攏，稱為桶狀畸變。

用數(shù)學(xué)公式來表示：

其中，X為相機中的坐標；X為真實世界坐標；K為內(nèi)參矩陣；RT為外參矩陣 K為內(nèi)參矩陣，是相機內(nèi)部參數(shù)組成的一個3*3的矩陣，其中，代表焦距；S為畸變參數(shù)為中心點坐標，a為縱橫比例參數(shù)，我們可以默認設(shè)為1，所以 RT為外參矩陣，R是描述照相機方向的旋轉(zhuǎn)矩陣，T是描述照相機中心位置的三維平移向量。

二、算法原理

1.整體流程

2.模型假設(shè)

3.模型求解

（1）內(nèi)外參數(shù)求解

我們令，則

其中，H為一個3*3的矩陣，并且有一個元素作為齊次坐標。因此，H有8個自由度。

現(xiàn)在有8個自由度需要求解，所以需要四個對應(yīng)點。也就是四個點就可以求出圖像平面到世界平面的單應(yīng)性矩陣H。

我想，張氏標定法選用的棋盤格作為標定板的原因除了角點方便檢測的另外一個原因可能就是這個吧。

通過4個點，我們就可以可以獲得單應(yīng)性矩陣H。但是，H是內(nèi)參陣和外參陣的合體。我們想要最終分別獲得內(nèi)參和外參。所以需要想個辦法，先把內(nèi)參求出來。然后外參也就隨之解出了。觀察一下這個式子：

我們可以知道以下約束條件:

①,R1R2 正交，也就是說 R1 R2=0。其實這個不難理解，因為 R1 R2 是分別繞x軸和y軸得到的，而x軸和y軸均垂直z軸。

②旋轉(zhuǎn)向量的模為1，也就是說R1=R2=1,這是因為旋轉(zhuǎn)不改變尺度。

根據(jù)這兩個約束條件，經(jīng)過數(shù)學(xué)變換，我們可以得到：

觀察上面的兩個式子，我們可以看出，由于H1和H2是通過單應(yīng)性求解出來的，所以我們要求解的參數(shù)就變成A矩陣中未知的5個參數(shù)。我們可以通過三個單應(yīng)性矩陣來求解這5個參數(shù)，利用三個單應(yīng)性矩陣在兩個約束下可以生成6個方程。其中，三個單應(yīng)性矩陣可以通過三張對同一標定板不同角度和高度的照片獲得。

用數(shù)學(xué)公式來表達如下：

我們很容易發(fā)現(xiàn)B是一個對稱陣，所以B的有效元素就剩下6個，即

進一步化簡：

通過計算，我們可以得到

利用上面提到的兩個約束條件，我們可以得到下面的方程組：

這個方程組的本質(zhì)和前面那兩個用h和A組成的約束條件方程組是一樣的。

通過至少含一個棋盤格的三幅圖像，應(yīng)用上述公式我們就可以估算出B了。得到B后，我們通過cholesky分解，就可以得到攝相機機的內(nèi)參陣A的六個自由度，即：

再根據(jù)化簡可得外部參數(shù)，即：

（2）畸變系數(shù)求解

在文章的開始，我們就講到真實的鏡頭并非理想的透視成像，而是帶有不同程度的畸變。理論上鏡頭的畸變包括徑向畸變和切向畸變，切向畸變影響較小，通常只考慮徑向畸變，而且在徑向畸變的求解中，僅考慮了起主導(dǎo)的二元泰勒級數(shù)展開的前兩個系數(shù)。

具體推導(dǎo)，參考https://blog.csdn.net/onthewaysuccess/article/details/40736177

（3）精度優(yōu)化

在張正友標定法中，使用了兩次極大似然估計策略，第一次是在不考慮畸變的情況下求解內(nèi)參和外參，第二次是求解實際的畸變系數(shù)。

極大似然參數(shù)估計，參考https://blog.csdn.net/onthewaysuccess/article/details/40717213

三、算法實現(xiàn)

1.main.py

#!usr/bin/env/ python

# _*_ coding:utf-8 _*_

import cv2 as cv

import numpy as np

import os

from step.homography import get_homography

from step.intrinsics import get_intrinsics_param

from step.extrinsics import get_extrinsics_param

from step.distortion import get_distortion

from step.refine_all import refinall_all_param

def calibrate():

#求單應(yīng)矩陣

H = get_homography(pic_points, real_points_x_y)

#求內(nèi)參

intrinsics_param = get_intrinsics_param(H)

#求對應(yīng)每幅圖外參

extrinsics_param = get_extrinsics_param(H, intrinsics_param)

#畸變矯正

k = get_distortion(intrinsics_param, extrinsics_param, pic_points, real_points_x_y)

#微調(diào)所有參數(shù)

[new_intrinsics_param, new_k, new_extrinsics_param] = refinall_all_param(intrinsics_param,

k, extrinsics_param, real_points, pic_points)

print("intrinsics_parm: ", new_intrinsics_param)

print("distortionk: ", new_k)

print("extrinsics_parm: ", new_extrinsics_param)

if __name__ == "__main__":

file_dir = r'..pic'

# 標定所用圖像

pic_name = os.listdir(file_dir)

# 由于棋盤為二維平面，設(shè)定世界坐標系在棋盤上，一個單位代表一個棋盤寬度，產(chǎn)生世界坐標系三維坐標

cross_corners = [9, 6] #棋盤方塊交界點排列

real_coor = np.zeros((cross_corners[0] * cross_corners[1], 3), np.float32)

real_coor[:, :2] = np.mgrid[0:9, 0:6].T.reshape(-1, 2)

real_points = []

real_points_x_y = []

pic_points = []

for pic in pic_name:

pic_path = os.path.join(file_dir, pic)

pic_data = cv.imread(pic_path)

# 尋找到棋盤角點

succ, pic_coor = cv.findChessboardCorners(pic_data, (cross_corners[0], cross_corners[1]), None)

if succ:

# 添加每幅圖的對應(yīng)3D-2D坐標

pic_coor = pic_coor.reshape(-1, 2)

pic_points.append(pic_coor)

real_points.append(real_coor)

real_points_x_y.append(real_coor[:, :2])

calibrate()

2.homography.py

這是用于求解單應(yīng)性矩陣的文件

#!usr/bin/env/ python

# _*_ coding:utf-8 _*_

import numpy as np

from scipy import optimize as opt

#求輸入數(shù)據(jù)的歸一化矩陣

def normalizing_input_data(coor_data):

x_avg = np.mean(coor_data[:, 0])

y_avg = np.mean(coor_data[:, 1])

sx = np.sqrt(2) / np.std(coor_data[:, 0])

sy = np.sqrt(2) / np.std(coor_data[:, 1])

norm_matrix = np.matrix([[sx, 0, -sx * x_avg],

[0, sy, -sy * y_avg],

[0, 0, 1]])

return norm_matrix

#求取初始估計的單應(yīng)矩陣

def get_initial_H(pic_coor, real_coor):

# 獲得歸一化矩陣

pic_norm_mat = normalizing_input_data(pic_coor)

real_norm_mat = normalizing_input_data(real_coor)

M = []

for i in range(len(pic_coor)):

#轉(zhuǎn)換為齊次坐標

single_pic_coor = np.array([pic_coor[i][0], pic_coor[i][1], 1])

single_real_coor = np.array([real_coor[i][0], real_coor[i][1], 1])

#坐標歸一化

pic_norm = np.dot(pic_norm_mat, single_pic_coor)

real_norm = np.dot(real_norm_mat, single_real_coor)

#構(gòu)造M矩陣

M.append(np.array([-real_norm.item(0), -real_norm.item(1), -1,

0, 0, 0,

pic_norm.item(0) * real_norm.item(0), pic_norm.item(0) * real_norm.item(1), pic_norm.item(0)]))

M.append(np.array([0, 0, 0,

-real_norm.item(0), -real_norm.item(1), -1,

pic_norm.item(1) * real_norm.item(0), pic_norm.item(1) * real_norm.item(1), pic_norm.item(1)]))

#利用SVD求解M * h = 0中h的解

U, S, VT = np.linalg.svd((np.array(M, dtype='float')).reshape((-1, 9)))

# 最小的奇異值對應(yīng)的奇異向量,S求出來按大小排列的，最后的最小

H = VT[-1].reshape((3, 3))

H = np.dot(np.dot(np.linalg.inv(pic_norm_mat), H), real_norm_mat)

H /= H[-1, -1]

return H

#返回估計坐標與真實坐標偏差

def value(H, pic_coor, real_coor):

Y = np.array([])

for i in range(len(real_coor)):

single_real_coor = np.array([real_coor[i, 0], real_coor[i, 1], 1])

U = np.dot(H.reshape(3, 3), single_real_coor)

U /= U[-1]

Y = np.append(Y, U[:2])

Y_NEW = (pic_coor.reshape(-1) - Y)

return Y_NEW

#返回對應(yīng)jacobian矩陣

def jacobian(H, pic_coor, real_coor):

J = []

for i in range(len(real_coor)):

sx = H[0]*real_coor[i][0] + H[1]*real_coor[i][1] +H[2]

sy = H[3]*real_coor[i][0] + H[4]*real_coor[i][1] +H[5]

w = H[6]*real_coor[i][0] + H[7]*real_coor[i][1] +H[8]

w2 = w * w

J.append(np.array([real_coor[i][0]/w, real_coor[i][1]/w, 1/w,

0, 0, 0,

-sx*real_coor[i][0]/w2, -sx*real_coor[i][1]/w2, -sx/w2]))

J.append(np.array([0, 0, 0,

real_coor[i][0]/w, real_coor[i][1]/w, 1/w,

-sy*real_coor[i][0]/w2, -sy*real_coor[i][1]/w2, -sy/w2]))

return np.array(J)

#利用Levenberg Marquart算法微調(diào)H

def refine_H(pic_coor, real_coor, initial_H):

initial_H = np.array(initial_H)

final_H = opt.leastsq(value,

initial_H,

Dfun=jacobian,

args=(pic_coor, real_coor))[0]

final_H /= np.array(final_H[-1])

return final_H

#返回微調(diào)后的H

def get_homography(pic_coor, real_coor):

refined_homographies =[]

error = []

for i in range(len(pic_coor)):

initial_H = get_initial_H(pic_coor[i], real_coor[i])

final_H = refine_H(pic_coor[i], real_coor[i], initial_H)

refined_homographies.append(final_H)

return np.array(refined_homographies)

3.intrinsics.py

這是用于求解內(nèi)參矩陣的文件

#!usr/bin/env/ python

# _*_ coding:utf-8 _*_

import numpy as np

#返回pq位置對應(yīng)的v向量

def create_v(p, q, H):

H = H.reshape(3, 3)

return np.array([

H[0, p] * H[0, q],

H[0, p] * H[1, q] + H[1, p] * H[0, q],

H[1, p] * H[1, q],

H[2, p] * H[0, q] + H[0, p] * H[2, q],

H[2, p] * H[1, q] + H[1, p] * H[2, q],

H[2, p] * H[2, q]

])

#返回相機內(nèi)參矩陣A

def get_intrinsics_param(H):

#構(gòu)建V矩陣

V = np.array([])

for i in range(len(H)):

V = np.append(V, np.array([create_v(0, 1, H[i]), create_v(0, 0 , H[i])- create_v(1, 1 , H[i])]))

#求解V*b = 0中的b

U, S, VT = np.linalg.svd((np.array(V, dtype='float')).reshape((-1, 6)))

#最小的奇異值對應(yīng)的奇異向量,S求出來按大小排列的，最后的最小

b = VT[-1]

#求取相機內(nèi)參

w = b[0] * b[2] * b[5] - b[1] * b[1] * b[5] - b[0] * b[4] * b[4] + 2 * b[1] * b[3] * b[4] - b[2] * b[3] * b[3]

d = b[0] * b[2] - b[1] * b[1]

alpha = np.sqrt(w / (d * b[0]))

beta = np.sqrt(w / d**2 * b[0])

gamma = np.sqrt(w / (d**2 * b[0])) * b[1]

uc = (b[1] * b[4] - b[2] * b[3]) / d

vc = (b[1] * b[3] - b[0] * b[4]) / d

return np.array([

[alpha, gamma, uc],

[0, beta, vc],

[0, 0, 1]

])

4.extrinsics.py

這是用于求解外參矩陣的文件

#!usr/bin/env/ python

# _*_ coding:utf-8 _*_

import numpy as np

#返回每一幅圖的外參矩陣[R|t]

def get_extrinsics_param(H, intrinsics_param):

extrinsics_param = []

inv_intrinsics_param = np.linalg.inv(intrinsics_param)

for i in range(len(H)):

h0 = (H[i].reshape(3, 3))[:, 0]

h1 = (H[i].reshape(3, 3))[:, 1]

h2 = (H[i].reshape(3, 3))[:, 2]

scale_factor = 1 / np.linalg.norm(np.dot(inv_intrinsics_param, h0))

r0 = scale_factor * np.dot(inv_intrinsics_param, h0)

r1 = scale_factor * np.dot(inv_intrinsics_param, h1)

t = scale_factor * np.dot(inv_intrinsics_param, h2)

r2 = np.cross(r0, r1)

R = np.array([r0, r1, r2, t]).transpose()

extrinsics_param.append(R)

return extrinsics_param

5.distortion.py

這是用于求解畸變矯正系數(shù)的文件

#!usr/bin/env/ python

# _*_ coding:utf-8 _*_

import numpy as np

#返回畸變矯正系數(shù)k0,k1

def get_distortion(intrinsic_param, extrinsic_param, pic_coor, real_coor):

D = []

d = []

for i in range(len(pic_coor)):

for j in range(len(pic_coor[i])):

#轉(zhuǎn)換為齊次坐標

single_coor = np.array([(real_coor[i])[j, 0], (real_coor[i])[j, 1], 0, 1])

#利用現(xiàn)有內(nèi)參及外參求出估計圖像坐標

u = np.dot(np.dot(intrinsic_param, extrinsic_param[i]), single_coor)

[u_estim, v_estim] = [u[0]/u[2], u[1]/u[2]]

coor_norm = np.dot(extrinsic_param[i], single_coor)

coor_norm /= coor_norm[-1]

#r = np.linalg.norm((real_coor[i])[j])

r = np.linalg.norm(coor_norm)

D.append(np.array([(u_estim - intrinsic_param[0, 2]) * r ** 2, (u_estim - intrinsic_param[0, 2]) * r ** 4]))

D.append(np.array([(v_estim - intrinsic_param[1, 2]) * r ** 2, (v_estim - intrinsic_param[1, 2]) * r ** 4]))

#求出估計坐標與真實坐標的殘差

d.append(pic_coor[i][j, 0] - u_estim)

d.append(pic_coor[i][j, 1] - v_estim)

'''

D.append(np.array([(pic_coor[i][j, 0] - intrinsic_param[0, 2]) * r ** 2, (pic_coor[i][j, 0] - intrinsic_param[0, 2]) * r ** 4]))

D.append(np.array([(pic_coor[i][j, 1] - intrinsic_param[1, 2]) * r ** 2, (pic_coor[i][j, 1] - intrinsic_param[1, 2]) * r ** 4]))

#求出估計坐標與真實坐標的殘差

d.append(u_estim - pic_coor[i][j, 0])

d.append(v_estim - pic_coor[i][j, 1])

'''

D = np.array(D)

temp = np.dot(np.linalg.inv(np.dot(D.T, D)), D.T)

k = np.dot(temp, d)

'''

#也可利用SVD求解D * k = d中的k

U, S, Vh=np.linalg.svd(D, full_matrices=False)

temp_S = np.array([[S[0], 0],

[0, S[1]]])

temp_res = np.dot(Vh.transpose(), np.linalg.inv(temp_S))

temp_res_res = np.dot(temp_res, U.transpose())

k = np.dot(temp_res_res, d)

'''

return k

6.refine_all.py

這是用于微調(diào)參數(shù)的文件

#!usr/bin/env/ python

# _*_ coding:utf-8 _*_

import numpy as np

import math

from scipy import optimize as opt

#微調(diào)所有參數(shù)

def refinall_all_param(A, k, W, real_coor, pic_coor):

#整合參數(shù)

P_init = compose_paramter_vector(A, k, W)

#復(fù)制一份真實坐標

X_double = np.zeros((2 * len(real_coor) * len(real_coor[0]), 3))

Y = np.zeros((2 * len(real_coor) * len(real_coor[0])))

M = len(real_coor)

N = len(real_coor[0])

for i in range(M):

for j in range(N):

X_double[(i * N + j) * 2] = (real_coor[i])[j]

X_double[(i * N + j) * 2 + 1] = (real_coor[i])[j]

Y[(i * N + j) * 2] = (pic_coor[i])[j, 0]

Y[(i * N + j) * 2 + 1] = (pic_coor[i])[j, 1]

#微調(diào)所有參數(shù)

P = opt.leastsq(value,

P_init,

args=(W, real_coor, pic_coor),

Dfun=jacobian)[0]

#raial_error表示利用標定后的參數(shù)計算得到的圖像坐標與真實圖像坐標點的平均像素距離

error = value(P, W, real_coor, pic_coor)

raial_error = [np.sqrt(error[2 * i]**2 + error[2 * i + 1]**2) for i in range(len(error) // 2)]

print("total max error: ", np.max(raial_error))

#返回拆解后參數(shù)，分別為內(nèi)參矩陣，畸變矯正系數(shù)，每幅圖對應(yīng)外參矩陣

return decompose_paramter_vector(P)

#把所有參數(shù)整合到一個數(shù)組內(nèi)

def compose_paramter_vector(A, k, W):

alpha = np.array([A[0, 0], A[1, 1], A[0, 1], A[0, 2], A[1, 2], k[0], k[1]])

P = alpha

for i in range(len(W)):

R, t = (W[i])[:, :3], (W[i])[:, 3]

#旋轉(zhuǎn)矩陣轉(zhuǎn)換為一維向量形式

zrou = to_rodrigues_vector(R)

w = np.append(zrou, t)

P = np.append(P, w)

return P

#分解參數(shù)集合，得到對應(yīng)的內(nèi)參，外參，畸變矯正系數(shù)

def decompose_paramter_vector(P):

[alpha, beta, gamma, uc, vc, k0, k1] = P[0:7]

A = np.array([[alpha, gamma, uc],

[0, beta, vc],

[0, 0, 1]])

k = np.array([k0, k1])

W = []

M = (len(P) - 7) // 6

for i in range(M):

m = 7 + 6 * i

zrou = P[m:m+3]

t = (P[m+3:m+6]).reshape(3, -1)

#將旋轉(zhuǎn)矩陣一維向量形式還原為矩陣形式

R = to_rotation_matrix(zrou)

#依次拼接每幅圖的外參

w = np.concatenate((R, t), axis=1)

W.append(w)

W = np.array(W)

return A, k, W

#返回從真實世界坐標映射的圖像坐標

def get_single_project_coor(A, W, k, coor):

single_coor = np.array([coor[0], coor[1], coor[2], 1])

#'''

coor_norm = np.dot(W, single_coor)

coor_norm /= coor_norm[-1]

#r = np.linalg.norm(coor)

r = np.linalg.norm(coor_norm)

uv = np.dot(np.dot(A, W), single_coor)

uv /= uv[-1]

#畸變

u0 = uv[0]

v0 = uv[1]

uc = A[0, 2]

vc = A[1, 2]

#u = (uc * r**2 * k[0] + uc * r**4 * k[1] - u0) / (r**2 * k[0] + r**4 * k[1] - 1)

#v = (vc * r**2 * k[0] + vc * r**4 * k[1] - v0) / (r**2 * k[0] + r**4 * k[1] - 1)

u = u0 + (u0 - uc) * r**2 * k[0] + (u0 - uc) * r**4 * k[1]

v = v0 + (v0 - vc) * r**2 * k[0] + (v0 - vc) * r**4 * k[1]

'''

uv = np.dot(W, single_coor)

uv /= uv[-1]

# 透鏡矯正

x0 = uv[0]

y0 = uv[1]

r = np.linalg.norm(np.array([x0, y0]))

k0 = 0

k1 = 0

x = x0 * (1 + r ** 2 * k0 + r ** 4 * k1)

y = y0 * (1 + r ** 2 * k0 + r ** 4 * k1)

#u = A[0, 0] * x + A[0, 2]

#v = A[1, 1] * y + A[1, 2]

[u, v, _] = np.dot(A, np.array([x, y, 1]))

'''

return np.array([u, v])

#返回所有點的真實世界坐標映射到的圖像坐標與真實圖像坐標的殘差

def value(P, org_W, X, Y_real):

M = (len(P) - 7) // 6

N = len(X[0])

A = np.array([

[P[0], P[2], P[3]],

[0, P[1], P[4]],

[0, 0, 1]

])

Y = np.array([])

for i in range(M):

m = 7 + 6 * i

#取出當前圖像對應(yīng)的外參

w = P[m:m + 6]

# 不用旋轉(zhuǎn)矩陣的變換是因為會有精度損失

'''

R = to_rotation_matrix(w[:3])

t = w[3:].reshape(3, 1)

W = np.concatenate((R, t), axis=1)

'''

W = org_W[i]

#計算每幅圖的坐標殘差

for j in range(N):

Y = np.append(Y, get_single_project_coor(A, W, np.array([P[5], P[6]]), (X[i])[j]))

error_Y = np.array(Y_real).reshape(-1) - Y

return error_Y

#計算對應(yīng)jacobian矩陣

def jacobian(P, WW, X, Y_real):

M = (len(P) - 7) // 6

N = len(X[0])

K = len(P)

A = np.array([

[P[0], P[2], P[3]],

[0, P[1], P[4]],

[0, 0, 1]

])

res = np.array([])

for i in range(M):

m = 7 + 6 * i

w = P[m:m + 6]

R = to_rotation_matrix(w[:3])

t = w[3:].reshape(3, 1)

W = np.concatenate((R, t), axis=1)

for j in range(N):

res = np.append(res, get_single_project_coor(A, W, np.array([P[5], P[6]]), (X[i])[j]))

#求得x, y方向?qū)[k]的偏導(dǎo)

J = np.zeros((K, 2 * M * N))

for k in range(K):

J[k] = np.gradient(res, P[k])

return J.T

#將旋轉(zhuǎn)矩陣分解為一個向量并返回，Rodrigues旋轉(zhuǎn)向量與矩陣的變換,最后計算坐標時并未用到，因為會有精度損失

def to_rodrigues_vector(R):

p = 0.5 * np.array([[R[2, 1] - R[1, 2]],

[R[0, 2] - R[2, 0]],

[R[1, 0] - R[0, 1]]])

c = 0.5 * (np.trace(R) - 1)

if np.linalg.norm(p) == 0:

if c == 1:

zrou = np.array([0, 0, 0])

elif c == -1:

R_plus = R + np.eye(3, dtype='float')

norm_array = np.array([np.linalg.norm(R_plus[:, 0]),

np.linalg.norm(R_plus[:, 1]),

np.linalg.norm(R_plus[:, 2])])

v = R_plus[:, np.where(norm_array == max(norm_array))]

u = v / np.linalg.norm(v)

if u[0] < 0 or (u[0] == 0 and u[1] < 0) or (u[0] == u[1] and u[0] == 0 and u[2] < 0):

u = -u

zrou = math.pi * u

else:

zrou = []

else:

u = p / np.linalg.norm(p)

theata = math.atan2(np.linalg.norm(p), c)

zrou = theata * u

return zrou

#把旋轉(zhuǎn)矩陣的一維向量形式還原為旋轉(zhuǎn)矩陣并返回

def to_rotation_matrix(zrou):

theta = np.linalg.norm(zrou)

zrou_prime = zrou / theta

W = np.array([[0, -zrou_prime[2], zrou_prime[1]],

[zrou_prime[2], 0, -zrou_prime[0]],

[-zrou_prime[1], zrou_prime[0], 0]])

R = np.eye(3, dtype='float') + W * math.sin(theta) + np.dot(W, W) * (1 - math.cos(theta))

return R

7.結(jié)果

拍攝的不同角度，不同高度的圖像

運行結(jié)果：

拍照的手機是華為p9，后置攝像頭是1200萬像素。

內(nèi)部參數(shù)矩陣是為：

[ 9.95397796e+02, -5.74043156e+00, 5.30659959e+02, 0.00000000e+00, 1.04963119e+03, 6.55565437e+02, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]

因為代碼是以一個方格為一個單位，沒有考慮單位長度，所以要求真實的參數(shù)應(yīng)該乘一個單位長度，博主采用的方格的尺寸是2.98cm的，自己拿excel畫的，get了一個新技能~~

審核編輯：郭婷