从标定图像估计一般相机矩阵

标定约束、归一化与相机可辨识性

第二章把相机矩阵当已知输入，本章用独立标定照片估它。标定靶是制造时已知三维位置的十六个彩色标记，坐标以米记录；未知珠子场景的坐标不提供给标定器。检测器从每张标定 PNG 自动提取区域，颜色描述子与标靶已知颜色匹配，才得到三维点与像素对。这样的已知几何是测量工具，不是偷用未知对象真值。标定靶误差也会进入结果，本合成器把它理想化为精确。

待估矩阵 P 有十二个元素，但整体非零比例不影响像素，因而一般投影相机有十一个独立比例自由度。每个三维到二维对应给两条独立线性约束，因此至少需要六个点，数量只是必要条件。若全体标定点共面，图像只能约束该平面到图像的单应关系，无法充分确定离平面的映射。程序检查点集减去均值后的秩为三，再进入 DLT；这不意味着所有接近共面的集合都数值稳定。

令三维齐次点为四维列向量，P 三行的转置分别记为列向量 p 一、p 二、p 三。透视式为 u 等于第一行点积除第三行点积，v 同理。把分母移到左边，得到第一行点积减 u 乘第三行点积为零，以及第二行点积减 v 乘第三行点积为零。把 P 的元素按行排成十二维未知向量，每个对应贡献两行设计矩阵 A。约束是齐次的，零矩阵虽满足方程，却没有相机意义。

在无噪声且配置充分时，设计矩阵有一维零空间；有检测误差时，通常没有严格零空间，改为寻找单位向量使代数残差最小。这里导入奇异值分解的数值结论：A 等于 U 乘非负对角奇异值再乘 V 转置，正交变换保持欧氏长度，因此最小残差对应最小奇异值的右奇异向量。取 V 转置最后一行并恢复三乘四形状。这个代数目标计算方便，却还不是原始像素重投影目标。

为什么先归一化？世界坐标可能以毫米写成数千，像素坐标可能几百，常数齐次项仍为一，设计矩阵各列尺度会差很远。有限精度线性代数的误差可能被放大，也会让代数目标以单位大小决定不同项权重。我们把二维点平移到均值零，再统一缩放，使平均欧氏半径为根号二；三维点使半径为根号三。这里是平均距离，不是距离平方平均后的根号，二者必须按实现一致。

反归一化顺序可以直接验证。将原世界齐次点先乘 T_X，归一化相机 P_n 将它送到归一化像素，再乘 T_x 的逆回到原像素。因此完整映射正是逆 T_x、P_n、T_X 的乘积。不能把两个变换都放到左侧。恢复后矩阵依然有整体尺度自由度；项目用第三行前三项的范数归一，再检查标靶的平均深度符号，必要时整体取负，固定前向约定。

完整数值构造：八个点怎样变成十六行方程

先把所有可复算数字固定下来。相机取P=[[2,0,0,0],[0,2,0,0],[0,0,1,0]]，像素主点为零、焦距为二，世界点单位为米。八个点是X、Y各取正负一，Z取三或五，故不是同一平面。这里只为手算采用小焦距；实际PNG项目焦距仍是95像素，不因这个算例而修改。

将每点齐次坐标q=[X,Y,Z,1]按行存储，A的u行是[q,0,0,0,0,−u q]，v行是[0,0,0,0,q,−v q]。例如第一个点的u行为[1,1,3,1,0,0,0,0,−2/3,−2/3,−2,−2/3]；v行为[0,0,0,0,1,1,3,1,−2/3,−2/3,−2,−2/3]。第二个点的u行为[−1,1,3,1,0,0,0,0,−2/3,2/3,2,2/3]，v行为[0,0,0,0,−1,1,3,1,2/3,−2/3,−2,−2/3]。这四行已把符号、行展开顺序和齐次常数明确到每个元素。

剩余六点逐行代入同一个规则，得到16×12矩阵。无噪声时它的秩应为11，零空间只有整体比例一个方向。P按行展开后的欧氏长度为三，因此SVD可能返回正或负的[2,0,0,0,0,2,0,0,0,0,1,0]/3。我们最后令第三行第三项为一，恢复表中的P，而不是把SVD向量第一项强行设为一。

这组点的三维均值为(0,0,4)，所有点距均值均为根号三，所以TX只将Z减四，不需额外缩放。像素均值为零，平均半径是8根号二/15，故Tx的两条像素轴缩放均为15/8。由Pn=Tx P TX逆，得到Pn=[[15/4,0,0,0],[0,15/4,0,0],[0,0,1,4]]，仍允许共同非零比例。恢复时必须左乘Tx逆、右乘TX。

import numpy as np
P=np.array([[2.,0,0,0],[0,2,0,0],[0,0,1,0]])
X=np.array([[x,y,z,1.] for z in [3.,5.] for y in [1.,-1.] for x in [1.,-1.]])
def pixels(camera,points):
    q=points@camera.T
    return q[:,:2]/q[:,2:3]
uv=pixels(P,X)
# 从正文的两行规则构造；不导入项目calibrate_dlt。
A=np.array([row for q,(u,v) in zip(X,uv)
            for row in [np.r_[q,np.zeros(4),-u*q],np.r_[np.zeros(4),q,-v*q]]])
assert np.linalg.matrix_rank(A)==11
_,singular,Vt=np.linalg.svd(A)
Q=Vt[-1].reshape(3,4);Q=Q/Q[2,2]
np.testing.assert_allclose(Q,P,atol=1e-12)
TX=np.eye(4);TX[2,3]=-4
Tx=np.diag([15/8,15/8,1.])
Pn=np.array([[15/4,0,0,0],[0,15/4,0,0],[0,0,1,4.]])
restored=np.linalg.inv(Tx)@Pn@TX
np.testing.assert_allclose(restored,P,atol=1e-14)
heldout=np.array([[.3,-.4,4.,1.]])
np.testing.assert_allclose(pixels(Q,heldout),[[.15,-.2]],atol=1e-12)
print('rank',np.linalg.matrix_rank(A),'smallest singular',singular[-1])
print('camera',Q.tolist(),'heldout',pixels(Q,heldout).tolist())

同平面伪解也给出两台具体矩阵

P′横向分子比P多Z/2−2。在Z=4时增量为零，因此平面内任意点都同投影；在X=1、Y=0、Z=5时，P给u=.4，P′给u=.5，差.1像素。两台相机都能完美解释同一平面观测，却对离平面点给不同答案，这才是“多收集同平面点仍无帮助”的具体自由度。

import numpy as np
P=np.array([[2.,0,0,0],[0,2,0,0],[0,0,1,0]])
Q=P.copy();Q[0,2]=.5;Q[0,3]=-2
plane=np.array([[1.,0,4,1],[-2.,1,4,1],[.3,-.7,4,1]])
def project(M,X):
    y=X@M.T
    return y[:,:2]/y[:,2:3]
np.testing.assert_allclose(project(P,plane),project(Q,plane),atol=1e-14)
other=np.array([[1.,0,5,1]])
assert abs(project(P,other)[0,0]-.4)<1e-14
assert abs(project(Q,other)[0,0]-.5)<1e-14
print('same-plane equality; heldout difference',project(Q,other)-project(P,other))

算例一：反归一化可逆性与单位变化

设二维归一化把像素减去一百、五十后乘零点零一，三维归一化把世界点减去中心后乘二。若已经求得归一化相机，把原三维点依次按这三个步骤映射，再除齐次尺度，应与直接用恢复的 P 相同。检验不需要知道相机真实姿态，只需抽任意未参加估计的点比较两个计算路径。如果结果不同，最先检查矩阵乘法顺序与行列存储，而不是重新调奇异值容差。

把全部世界标靶坐标从米改写为毫米，相机矩阵前三列应相应除以一千，平移列相对保持，使齐次乘积不变；整体归一后的具体矩阵数值还可再乘共同因子。重投影不变并不意味着三角化输出的单位不变，新输出将是毫米。项目把输入单位和坐标约定一起保存，之后比较三维误差必须先统一单位，不能把米与毫米的差当算法精度提高一千倍。

估出 P 后要计算几何残差：将每个已知三维点重投影到像素，再减检测位置，二维长度平方平均后开根号。这个 RMS 使用每点二维欧氏误差，不能与每坐标 MSE 的根号混同。还要保存逐点误差与图像位置，检查残差是否在边缘较大、是否有少数错误匹配被平均遮住。项目五张标定图的 RMS 约零点零二二到零点零二九像素，是受控干净标记的实测，不是现实相机精度承诺。

算例二：代数残差小仍能有错误相机

假设所有已知点深度都为四米。齐次三维点的第三分量总是第四分量的四倍，设计矩阵相关列无法区分某些相机系数。可以改变 P 的深度列并相应补偿平移列，让平面上的像素完全不变，却让深度五米的未知点投影大幅改变。这不是多收集几千个同平面点就能解决的独立信息缺失。需要多姿态平面标定加内参约束，或本项目采用的非共面靶来恢复一般 P。

另一个独立退化来自像素侧：世界点虽非共面，若错误观测全部落在图像同一条直线上，就可能得到秩二映射，而非合法的秩三透视相机。新增实现先要求N×3与N×2行数一致、所有数字有限，再检查归一化世界点和二维像素的空间展开；设计矩阵的次小奇异值必须与最大值保持可分辨比例，最后检查恢复相机左3×3块的相对奇异值。不能只看最小奇异值接近零，因为合法齐次解本来就要有一个零方向，退化的危险是出现第二个未约束方向。

本教学实现明确使用归一化空间比例阈值：世界/像素展开和设计次小奇异值为1e−8，相机左块为1e−10。它们是数值可辨识性策略，不是所有噪声模型下的统计最优阈值。归一化减少米/毫米和像素尺度影响；直接判断未归一化行列式“小于某常数”会把单位变化误当退化。新增回归包含非共面世界加共线像素、重复观测、近共面/近共线、错行数与NaN，并确认原正常标定相机不变。

即便不严格共面，若标定点只在非常薄的深度层里变化，设计矩阵也可能病态。末尾两个奇异值的相对大小、标靶空间覆盖、对像素扰动的敏感性都应一起看。最小奇异值很小只说明有一个近似解，并不保证该解对噪声稳定。固定种子给检测位置加零点一像素扰动，重新标定，再看留出点重投影变化，是一种具体压力测试；它给数值证据，不是普适误差界。

一般 P 估计不直接保证其左三乘三部分具有期望的零斜切内参形式。若要提取 K、R、t，可对左块做 RQ 分解，调整对角符号和旋转行列式，再设内参末元素为一；这是另一条有约束的相机解释步骤。本项目的三角化只需要 P，受控渲染另使用公开的焦距尺度，并未实现真实镜头的全内外参与畸变联合非线性标定。不能把 DLT 输出的任意矩阵默认为所有现实约束均满足。

进一步精化通常最小化像素重投影平方和，而不是代数残差。把相机用合法旋转、平移和内参参数化，对参数求导后迭代，才不会在优化中任意离开旋转流形。非线性目标可能有局部极小值，所以 DLT 的作用是给合理初值，不是保证全局最优。本课程将该扩展的接口与原理讲清，当前项目不冒称运行了 OpenCV 标定器或全场景束调整。

程序会对五个相机分别读取各自标定图片，包括保留场景视角的标定图。保留视角的相机参数来自独立标靶，并不使用保留场景图像来调整珠子位置。这对应“已标定新视角”的评价条件，而不是完全未知相机姿态的重建。评价时必须把这个条件写出来，否则读者可能以为项目还实现了任意视角定位。未知姿态可以作为后续求解问题，但不在当前实测范围。

源数据的可审计性与数值稳定性是互补问题。清单记录每个 PNG 和标靶文件的散列，能够确认运行读的是什么；秩、奇异值和残差能够检查这些输入是否足以支持估计。散列正确的错误对应仍会给错误相机，残差很小的结果也可能来自错误单位。验收必须沿着文件、观测、对应、矩阵、残差和单位逐层检查，不让某一个通过标志遮蔽其他层的失败。

import numpy as np
from vision3d.core import look_at,project,calibrate_dlt
rng=np.random.default_rng(23)
X=rng.uniform([-1,-1,3],[1,1,5],(20,3))
P,_=look_at([.8,-.2,.1])
u,_=project(X,P)
Q,s=calibrate_dlt(X,u)
np.testing.assert_allclose(project(X,Q)[0],u,atol=1e-10)
print('reprojection max',abs(project(X,Q)[0]-u).max())

# design_rows
def design_rows(world,pixels):
    raise NotImplementedError('dlt')

# calibrate
def calibrate(world,pixels):
    raise NotImplementedError('dlt')

# design_rows
def design_rows(world,pixels):
 X=np.c_[np.asarray(world,dtype=float),np.ones(len(world))];uv=np.asarray(pixels,dtype=float);rows=[]
 for point,(u,v)in zip(X,uv):rows.extend([np.r_[point,np.zeros(4),-u*point],np.r_[np.zeros(4),point,-v*point]])
 return np.asarray(rows)

# calibrate
def calibrate(world,pixels):
 world=np.asarray(world,dtype=float);pixels=np.asarray(pixels,dtype=float)
 if world.shape[0]<6 or world.shape[1:]!=(3,)or pixels.shape!=(len(world),2)or not np.isfinite(world).all()or not np.isfinite(pixels).all():raise ValueError('finite Nx3/Nx2, N>=6')
 if np.linalg.matrix_rank(np.c_[world,np.ones(len(world))])<4:raise ValueError('coplanar input')
 A=design_rows(world,pixels);u,s,v=np.linalg.svd(A,full_matrices=False)
 if np.linalg.matrix_rank(A,tol=s[0]*1e-10)<11:raise ValueError('unidentifiable camera')
 P=v[-1].reshape(3,4)
 if np.linalg.matrix_rank(P,tol=np.linalg.norm(P)*1e-10)<3:raise ValueError('rank deficient camera')
 return P/np.linalg.norm(P)