单视图几何标定

参考链接：https://www.cv-xueba.club/

非畸变情况的单视图标定

• 求解$PM=0$的最小二乘解，由于是超定齐次方程组，通过SVD求解得到M
• 将M分解为$A[a_1,a_2,a_3]$和b，先计算内参数得到K矩阵，再计算R和T外参数矩阵。
  注意：最少选取六对点作为输入参数，这六对点不能位于同一平面因为三个点可以确定一个平面，第四个同平面的点可以被三个点表示，计算结果退化。

把图中黑白块单元格边长看作一个像素，手动输入计算参数。

Python程序

运行结果可以在colab中观看：colab:https://colab.research.google.com/drive/1HQN0BF0VmRP79pboOk3peWTqm6GewpaJ?usp=sharing

完全参考自 https://github.com/CV-xueba/Total3DExercises/tree/main/MVGlab01_camera-calibration-master

import numpy as np
from numpy import linalg as LA
class SingleCamera:
  # 初始化参数
  def __init__(self, world_coor, pixel_coor, n):
    self.__world_coor = world_coor
    self.__pixel_coor = pixel_coor
    self.__point_num = n

    self.__P = np.empty([2*self.__point_num,12],dtype=float)
    self.__roM = np.empty([3,4], dtype=float)
    self.__A = np.empty([3,3], dtype=float)
    self.__b = np.empty([3,1], dtype=float)
    self.__K = np.empty([3,3], dtype=float)
    self.__R = np.empty([3,3], dtype=float)
    self.__T = np.empty([3,3], dtype=float)
    
  # 构造P矩阵的形式 2n*12
  def composeP(self):
    P = np.empty([2*self.__point_num,12],dtype=float)
    i=0
    while(i< self.__point_num):
      Pi = self.__world_coor[i]
      fourzero =np.zeros([4],dtype=float)
      uminusPi = -self.__pixel_coor[i][0]*Pi
      vminusPi = -self.__pixel_coor[i][1]*Pi
      P[2*i] = np.hstack((Pi,fourzero,uminusPi))
      P[2*i+1] = np.hstack((fourzero,Pi,vminusPi))  
      i+=1
    print("Now P is with the form of :")
    print(P)
    print('\n')
    self.__P = P
    
  # 对PM=0中P奇异值分解得到M[A，b]
  def svdP(self):
    U,sigma,Vt = LA.svd(self.__P)
    V = np.transpose(Vt)
    # print(V[:,-1])
    roM = V[:,-1].reshape(3,4)
    print("roM:")
    print(roM)
    print('\n')
    A = roM[:,0:3].copy()
    b = roM[:,3:4].copy()
    print("A(3x3):")
    print(A)
    print("b(3x1):")
    print(b)
    self.__A = A
    self.__b = b
    self.__roM = roM
  
  # 计算内参数/外参数
  def solve(self):
    a3T = self.__A[2]
    a3Norm = LA.norm(a3T)
    rho = 1/a3Norm
    print("The rho is %f \n" % rho)

    a1T = self.__A[0]
    a2T = self.__A[1]
    cx = rho*rho*(np.dot(a1T,a3T))
    cy = rho*rho*(np.dot(a2T,a3T))
    print("The cx is %f \n The cy is %f" %(cx,cy))

    aCross13=np.cross(a1T,a3T)
    aCross23=np.cross(a2T,a3T)
    theta =np.arccos( -np.dot(aCross13,aCross23)/LA.norm(aCross13)/LA.norm(aCross23))
    print("The theta is %f: \n" %theta)

    alpha = rho*rho*LA.norm(aCross13)*np.sin(theta)
    beta = rho*rho*LA.norm(aCross23)*np.sin(theta)
    print("alpha:%f \n beta:%f \n"%(alpha,beta))

    # the K 
    K = np.array([alpha,-alpha*np.cos(theta),cx,0,beta/np.sin(theta),cy,0,0,1])
    K = K.reshape(3,3)
    print("The K:")
    print(K)
    print('\n')
    self.__K = K

    # compute R
    r1 = aCross23/LA.norm(aCross23)
    r3 = rho*a3T
    r2 = np.cross(r1,r3)
    R = np.hstack((r1,r2,r3))
    R.reshape(3,3)
    print("The Rotation:");
    print(R)
    print('\n')
    self.__R=R

    # compute T
    T = rho*np.dot(LA.inv(K),self.__b)
    print("The T:")
    print(T)
    print('\n')
    self.__T=T

  # 验证参数矩阵效果
  def selfcheck(self,w_check,c_check):
    my_size = c_check.shape[0]
    my_err = np.empty([my_size])

    for i in range(my_size):
      test_pix = np.dot(self.__roM,w_check[i])
      u = test_pix[0]/test_pix[2]
      v = test_pix[1]/test_pix[2]
      u_c = c_check[i][0]
      v_c = c_check[i][1]
      print("test point index:%d \n"%i)
      print("the predcit u & v:(%f,%f)\n"%(u,v))
      print("the orign u & v in world\n:(%f,%f)"%(u_c,v_c))
      my_err[i]=(abs(u-u_c)/u_c+abs(v-v_c)/v_c)/2
    average_err = my_err.sum()/my_size
    print("The average error is %f:\n"%average_err)

    if average_err > 0.1:
      print("average error is more than 0.1, is bad")
    else:
      print("average error is smaller tha 0.1, is acceptable")

测试结果:可以看到计算结果相对准确。

test point index:0 

the predcit u & v:(369.955752,296.987629)

the orign u & v in world
:(369.000000,297.000000)
test point index:1 

the predcit u & v:(530.973337,482.573653)

the orign u & v in world
:(531.000000,484.000000)
test point index:2 

the predcit u & v:(639.757701,468.574861)

the orign u & v in world
:(640.000000,468.000000)
test point index:3 

the predcit u & v:(647.491900,332.448606)

the orign u & v in world
:(646.000000,333.000000)
test point index:4 

the predcit u & v:(559.840202,194.248688)

the orign u & v in world
:(556.000000,194.000000)
The average error is 0.001939:

average error is smaller tha 0.1, is acceptable

C++程序