如何不使用cv2.calibrateCamera实现相机标定？

相机标定的核心是求解内参矩阵K、畸变系数D，以及每张标定板对应的外参（旋转向量r_vec、平移向量t_vec），核心基于针孔相机模型与畸变模型。以下是手动实现的完整流程：

核心步骤与代码实现 #

1. 准备基础数据 #

和使用cv2.calibrateCamera时一致，需要：

• 世界坐标系下的3D点集合（pts3d_list，每个元素为单张图的N×3数组）
• 对应图像中的2D点集合（pts2d_list，每个元素为单张图的N×2数组）
• 图像尺寸（img_size，即grayColor.shape[::-1]）

2. 实现单应性矩阵求解（DLT算法） #

单应性矩阵描述3D平面到2D图像的映射关系，是求解内参的基础：

import numpy as np

def compute_homography(pts3d, pts2d):
    # 转换为齐次坐标
    pts3d_h = np.hstack([pts3d, np.ones((len(pts3d), 1))])
    pts2d_h = np.hstack([pts2d, np.ones((len(pts2d), 1))])
    
    A = []
    for i in range(len(pts3d)):
        X, Y, Z, _ = pts3d_h[i]
        u, v, _ = pts2d_h[i]
        A.append([-X, -Y, -Z, -1, 0, 0, 0, 0, u*X, u*Y, u*Z, u])
        A.append([0, 0, 0, 0, -X, -Y, -Z, -1, v*X, v*Y, v*Z, v])
    A = np.array(A)
    
    # SVD求解单应性矩阵
    U, S, Vt = np.linalg.svd(A)
    H = Vt[-1].reshape(3,4)
    # 归一化，保证H[2,3]=1
    H /= H[2,3]
    return H

3. 从单应性矩阵求解内参初始值（张正友标定法） #

通过多张图的单应性矩阵，构建内参的约束方程，求解初始内参K：

def build_intrinsic_constraints(H_list):
    constraints = []
    for H in H_list:
        h1 = H[:,0].reshape(3,1)
        h2 = H[:,1].reshape(3,1)
        # 构建内参约束的向量形式
        v12 = np.array([h1[0]*h2[0], h1[0]*h2[1]+h1[1]*h2[0], h1[1]*h2[1], 
                        h1[0]*h2[2]+h1[2]*h2[0], h1[1]*h2[2]+h1[2]*h2[1], h1[2]*h2[2]])
        constraints.append(v12)
        
        v11_v22 = np.array([h1[0]**2 - h2[0]**2, 2*(h1[0]*h1[1] - h2[0]*h2[1]), 
                            h1[1]**2 - h2[1]**2, 2*(h1[0]*h1[2] - h2[0]*h2[2]), 
                            2*(h1[1]*h1[2] - h2[1]*h2[2]), h1[2]**2 - h2[2]**2])
        constraints.append(v11_v22)
    
    constraints = np.array(constraints)
    # SVD求解内参的中间矩阵B
    U, S, Vt = np.linalg.svd(constraints)
    b = Vt[-1]
    B11, B12, B22, B13, B23, B33 = b
    
    # 从B矩阵恢复内参K
    v0 = (B12*B13 - B11*B23)/(B11*B22 - B12**2)
    lambd = B33 - (B13**2 + v0*(B12*B13 - B11*B23))/B11
    alpha = np.sqrt(lambd/B11)
    beta = np.sqrt(lambd*B11/(B11*B22 - B12**2))
    gamma = -B12*alpha**2*beta/lambd
    u0 = gamma*v0/alpha - B13*alpha**2/lambd
    
    K = np.array([[alpha, gamma, u0],
                  [0, beta, v0],
                  [0, 0, 1]])
    return K

4. 分解单应性矩阵得到外参初始值 #

利用已求得的内参K，从单应性矩阵中分解出旋转矩阵与平移向量，并转换为旋转向量：

def decompose_homography(H, K):
    K_inv = np.linalg.inv(K)
    h1 = H[:,0]
    h2 = H[:,1]
    h3 = H[:,2]
    
    # 计算尺度因子
    lambd = 1/np.linalg.norm(K_inv @ h1)
    r1 = lambd * K_inv @ h1
    r2 = lambd * K_inv @ h2
    r3 = np.cross(r1, r2)
    t = lambd * K_inv @ h3
    
    # 修正旋转矩阵的正交性（保证行列式为1）
    U, S, Vt = np.linalg.svd(np.column_stack([r1, r2, r3]))
    R = U @ Vt
    if np.linalg.det(R) < 0:
        R[:,2] *= -1
        t *= -1
    
    return R, t

5. 非线性优化（最小化重投影误差） #

通过Levenberg-Marquardt算法，优化内参、畸变系数与外参，最小化3D点到2D点的重投影误差：

from scipy.optimize import leastsq

def reprojection_error(params, K, pts3d_list, pts2d_list):
    k1, k2, p1, p2, k3 = params[:5]
    extrinsics = params[5:].reshape(-1,6)
    
    total_error = []
    for idx in range(len(pts3d_list)):
        pts3d = pts3d_list[idx]
        pts2d = pts2d_list[idx]
        # 提取外参并转换为旋转矩阵
        r_vec = extrinsics[idx,:3].reshape(3,1)
        t_vec = extrinsics[idx,3:].reshape(3,1)
        R, _ = cv2.Rodrigues(r_vec)
        
        # 投影到归一化图像坐标
        pts3d_h = np.hstack([pts3d, np.ones((len(pts3d),1))])
        proj_pts = (R @ pts3d_h[:,:3].T + t_vec).T
        x = proj_pts[:,0]/proj_pts[:,2]
        y = proj_pts[:,1]/proj_pts[:,2]
        
        # 应用畸变模型
        r2 = x**2 + y**2
        x_distorted = x*(1 + k1*r2 + k2*r2**2 + k3*r2**3) + 2*p1*x*y + p2*(r2 + 2*x**2)
        y_distorted = y*(1 + k1*r2 + k2*r2**2 + k3*r2**3) + p1*(r2 + 2*y**2) + 2*p2*x*y
        
        # 投影到图像像素坐标
        u = K[0,0]*x_distorted + K[0,1]*y_distorted + K[0,2]
        v = K[1,1]*y_distorted + K[1,2]
        
        # 记录误差
        total_error.extend(np.hstack([u - pts2d[:,0], v - pts2d[:,1]]))
    
    return np.array(total_error)

def non_linear_optimization(K, pts3d_list, pts2d_list, r_vecs_init, t_vecs_init):
    # 初始化优化参数：畸变系数初始为0，外参用初始值
    k_init = np.zeros(5)
    extrinsics_init = []
    for r_vec, t_vec in zip(r_vecs_init, t_vecs_init):
        extrinsics_init.extend(r_vec.flatten())
        extrinsics_init.extend(t_vec.flatten())
    extrinsics_init = np.array(extrinsics_init)
    params_init = np.hstack([k_init, extrinsics_init])
    
    # 执行优化
    params_opt, _ = leastsq(reprojection_error, params_init, args=(K, pts3d_list, pts2d_list))
    
    # 提取优化结果
    k1, k2, p1, p2, k3 = params_opt[:5]
    distortion = np.array([k1, k2, p1, p2, k3])
    extrinsics_opt = params_opt[5:].reshape(-1,6)
    r_vecs_opt = []
    t_vecs_opt = []
    for extrinsic in extrinsics_opt:
        r_vecs_opt.append(extrinsic[:3].reshape(3,1))
        t_vecs_opt.append(extrinsic[3:].reshape(3,1))
    
    return K, distortion, r_vecs_opt, t_vecs_opt

6. 整合完整标定流程 #

def custom_calibrate(pts3d_list, pts2d_list, img_size):
    # 计算所有单应性矩阵
    H_list = []
    for pts3d, pts2d in zip(pts3d_list, pts2d_list):
        H = compute_homography(pts3d, pts2d)
        H_list.append(H)
    
    # 求解内参初始值
    K_init = build_intrinsic_constraints(H_list)
    
    # 求解外参初始值
    r_vecs_init = []
    t_vecs_init = []
    for H in H_list:
        R, t = decompose_homography(H, K_init)
        r_vec, _ = cv2.Rodrigues(R)
        r_vecs_init.append(r_vec)
        t_vecs_init.append(t)
    
    # 非线性优化
    K_opt, distortion_opt, r_vecs_opt, t_vecs_opt = non_linear_optimization(K_init, pts3d_list, pts2d_list, r_vecs_init, t_vecs_init)
    
    # 计算平均重投影误差作为返回的ret值
    ret = np.mean(np.abs(reprojection_error(np.hstack([distortion_opt, np.array(r_vecs_opt).flatten(), np.array(t_vecs_opt).flatten()]), K_opt, pts3d_list, pts2d_list)))
    
    return ret, K_opt, distortion_opt, r_vecs_opt, t_vecs_opt

注意事项 #

• 上述代码依赖numpy、scipy，以及cv2.Rodrigues（若完全不想依赖OpenCV，可自行实现旋转向量与矩阵的转换）。
• 张正友标定法要求至少3张标定图，实际建议使用10-20张以提升标定精度。
• 若需要更高精度，可修改非线性优化部分，将内参的元素（alpha、beta等）加入优化变量。

内容的提问来源于stack exchange，提问作者Sina