【多视几何】两视图单应变换推导

变换层次

常见2d变换类型包含

• 欧式变换（4dof）旋转+平移

• 相似变换（5dof）旋转+平移+全局缩放

• 仿射变换（6dof）旋转+平移+两垂直方向非均匀缩放
  共线/平行线的长度比不变

• 射影/单应/投影变换 8dof
  物平面在像平面成像 即为一种单应
  平面引导的两视图间的映射关系 可以看作是单应的复合
  单应为数学意义上的群，因此多个单应复合后仍为单应

  射影不变量 直线上交比（cross ratio）即长度比之比不变
  共线且不同的四点 a,b,c,d 定义的交比为|ab|·|cd|/|ac|/|bd|

  射影变换和仿射变换最大的不同在于3×3变换矩阵第三行的前两个元素非0，这样导致射影空间下的理想点/无穷远点(x1,x2,0)可以映射到有限点。

物平面引导的两视单应

• 问题陈述
  已知两张影像各自的内外参数矩阵 K1,R1,C1 和 K2,R2,C2，则世界坐标系下一点P在像面成像坐标分别为 x = K1 · (R1 · P - C1) 和x = K2 · (R1 · P - C1)。已知世界坐标系下的某参数化平面n.tr · P + d = 0 其中n为单位法向

  如果是已知相机坐标系下的参数n, d 参考这篇博客 或转到世界坐标系下进行后续步骤

  求解两张影像的单应变换H（其中x2 = H · x1 表示）

• 公式推导
  影像1上二维像点 x1 （加齐次项1）生成无尺度三维物点 P = R1.inv · K1.inv · x1 + C1（差尺度因子）

  x_cam = K1.inv · x1
P = R1.inv · x_cam + C1

  重投影到影像2，生成二维像点 x2 = K2 · R2 · (P - C2)（需要将齐次项缩放至1）
  尺度不确定的情况下，两像之间的关系由基本矩阵F确定，为点到线的映射
  若三维物点在落在物方平面(n,d)上，则尺度确定，可以计算出点到点的映射关系，即单应矩阵H

  平面约束 n.tr · (R1.inv · K1.inv · x1 + C1) + d = 0
  分离常数项 n.tr · R1.inv · K1.inv · x1 + n.tr · C1 + d = 0
  两边同除常数项 n.tr / (n.tr · C1 + d) · R1.inv · K1.inv · x1 = -1
  最终简化为 常数 × 矩阵 × 变量 = 1

  两像变换 x2 = H · x1
  用参数矩阵表示 x2 = K2 · R2 · ( R1.inv · K1.inv · x1 + C1 - C2 )
  利用1代换x2 = K2 · R2 · ( R1.inv ·K1.inv · x1 - (C1 - C2) · n.tr / (n.tr · C1 + d) · R1.inv · K1.inv · x1 )
  分离参数 x2 = K2 · R2 · ( R1.inv · K1.inv - (C1 - C2) · n.tr / (n.tr · C1 + d) · R1.inv · K1.inv · x1
  合并同类项 x2 = K2 · R2 · ( Identity - (C1 - C2) · n.tr / (n.tr · C1 + d) ) · R1.inv · K1.inv· x1

  最终算得单应矩阵H = K2·R2·(Identity - (C1 - C2)·n.tr/(n.tr· C1 + d))·R1.inv·K1.inv·x1

• 代码

  cv::Matx33d FaceData::homograghyFromPairImage(Image& im1, Image& im2) {
  	const cv::Matx33d& K1 = im1.camera.K;
  	const cv::Matx33d& R1 = im1.camera.R;
  	const cv::Matx31d& C1 = im1.camera.C;
  
  	const cv::Matx33d& K2 = im2.camera.K;
  	const cv::Matx33d& R2 = im2.camera.R;
  	const cv::Matx31d& C2 = im2.camera.C;
  
  	const cv::Matx31d& normal(faceData.normal);  // 世界坐标系下的法向
  	const double d(faceData.d);                  // n.t() * P + d = 0
  
  	double param(n.t() * C1 + d);
   	cv::Matx31d t(C1 - C2); 
  	cv::Matx33d H(R2 * (cv::Matx33d::eye() - t * n.t() * INVERT(param)) * R1.t());
  	return K2 * H * K1.inv();
  
  void main() { 
  	// FaceData face; Image im1,im2;
  	cv::Matx33d H = face.homographyFromPairImage(im1, im2);
  	// x1 from im1, x2 is the correspondence in im2
  	x2 = H * x1;
  }