为您推荐:
精华内容
最热下载
问答
  • 5星
    3KB qq_49747343 2020-12-04 09:11:20
  • 5星
    12.28MB qq_32688843 2021-06-07 17:54:37
  • 5星
    7KB tkl32172 2021-04-19 20:04:50
  • 5星
    17.04MB weixin_43557432 2021-03-14 10:29:07
  • 4星
    39.09MB haojunjun 2021-08-01 19:29:38
  • 5星
    6.21MB GZXGYZ 2021-03-18 16:11:34
  • 5星
    223KB m0_52957036 2020-03-15 16:14:37
  • 5星
    15MB GZXGYZ 2021-03-18 14:46:04
  • 5星
    8.67MB GISuuser 2021-09-17 19:32:38
  • 5星
    53.55MB jiangwenping2009 2021-07-29 21:17:08
  • 4星
    1018KB weixin_42696333 2021-09-10 21:08:27
  • 37.1MB zhanzhanzhanbuaa 2018-06-11 18:55:14
  • 1.68MB hellodazhi 2021-09-30 00:00:20
  • 489KB weixin_38666785 2020-10-16 08:51:21
  • 4星
    56.03MB yanpeng_love 2018-08-14 12:45:37
  • OpenCV+OpenGL 双目立体视觉三维重建 OpenCV 双目测距(双目标定、双目校正和立体匹配) 真实场景的双目立体匹配(Stereo Matching)获取深度图详解 图像获取 双目相机拍摄获取 左右目图像 双目标定 ...

    Overview

    双目立体视觉的整体流程包括:

    • 图像采集
    • 双目标定
    • 双目矫正
    • 立体匹配
    • 三维重建

    1. 图像采集

    双目相机采集 左右目图像

    2. 双目标定

    通过 双目标定工具 对双目相机进行标定,得到如下结果参数:

    内参外参
    相机矩阵 K 1 , K 2 K_1, K_2 K1,K2旋转矩阵 R R R
    畸变系数 D 1 , D 2 D_1, D_2 D1,D2平移向量 t t t

    《Learning OpenCV》中对于 Translation 和 Rotation 的图示是这样的:

    示例代码:

    cv::Matx33d K1, K2, R;
    cv::Vec3d T;
    cv::Vec4d D1, D2;
    
    int flag = 0;
    flag |= cv::fisheye::CALIB_RECOMPUTE_EXTRINSIC;
    flag |= cv::fisheye::CALIB_CHECK_COND;
    flag |= cv::fisheye::CALIB_FIX_SKEW;
    
    cv::fisheye::stereoCalibrate(
            obj_points_, img_points_l_, img_points_r_,
            K1, D1, K2, D2, img_size_, R, T,
            flag, cv::TermCriteria(3, 12, 0));
    

    3. 双目矫正

    双目矫正 主要包括两方面:畸变矫正立体矫正

    利用 OpenCV的函数,主要分为

    • stereoRectify
    • initUndistortRectifyMap
    • remap

    stereoRectify

    根据双目标定的结果 K 1 , K 2 , D 1 , D 2 , R , t K_1, K_2, D_1, D_2, R, t K1,K2,D1,D2,R,t,利用 OpenCV函数 stereoRectify,计算得到如下参数

    • 左目 矫正矩阵(旋转矩阵) R 1 R_1 R1 (3x3)
    • 右目 矫正矩阵(旋转矩阵) R 2 R_2 R2 (3x3)
    • 左目 投影矩阵 P 1 P_1 P1 (3x4)
    • 右目 投影矩阵 P 2 P_2 P2 (3x4)
    • disparity-to-depth 映射矩阵 Q Q Q (4x4)

    其中,

    左右目投影矩阵(horizontal stereo, c x 1 ′ = c x 2 ′ {c_x}_1'={c_x}_2' cx1=cx2 if CV_CALIB_ZERO_DISPARITY is set)

    P 1 = [ f ′ 0 c x 1 ′ 0 0 f ′ c y ′ 0 0 0 1 0 ] P_1 = \begin{bmatrix} f' & 0 & {c_x}_1' & 0 \\ 0 & f' & c_y' & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} P1=f000f0cx1cy1000

    P 2 = [ f ′ 0 c x 2 ′ t x ′ ⋅ f ′ 0 f ′ c y ′ 0 0 0 1 0 ] P_2 = \begin{bmatrix} f' & 0 & {c_x}_2' & t_x' \cdot f' \\ 0 & f' & c_y' & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} P2=f000f0cx2cy1txf00

    where

    t x ′ = − B t_x' = -B tx=B

    disparity-to-depth 映射矩阵

    Q = [ 1 0 0 − c x 1 ′ 0 1 0 − c y ′ 0 0 0 f ′ 0 0 − 1 t x ′ c x 1 ′ − c x 2 ′ t x ′ ] Q = \begin{bmatrix} 1 & 0 & 0 & -{c_x}_1' \\ 0 & 1 & 0 & -c_y' \\ 0 & 0 & 0 & f' \\ 0 & 0 & -\frac{1}{t_x'} & \frac{ {c_x}_1'-{c_x}_2'}{t_x'} \end{bmatrix} Q=10000100000tx1cx1cyftxcx1cx2

    通过 P 2 P_2 P2 可计算出 基线 长度:

    baseline = B = − t x ′ = − P 2 ( 03 ) f ′ \begin{aligned} \text{baseline} = B = - t_x' = - \frac{ {P_2}_{(03)} }{f'} \end{aligned} baseline=B=tx=fP2(03)

    示例代码:

    cv::Mat R1, R2, P1, P2, Q;
    cv::fisheye::stereoRectify(
            K1, D1, K2, D2, img_size_, R, T,
            R1, R2, P1, P2, Q,
            CV_CALIB_ZERO_DISPARITY, img_size_, 0.0, 1.1);
    

    CameraInfo DKRP

    参考:sensor_msgs/CameraInfo Message

    • D: distortion parameters.

      • For “plumb_bob”, the 5 parameters are: (k1, k2, t1, t2, k3)
    • K: Intrinsic camera matrix for the raw (distorted) images.

      • Projects 3D points in the camera coordinate frame to 2D pixel coordinates using the focal lengths (fx, fy) and principal point (cx, cy).
        K = [ f x 0 c x 0 f y c y 0 0 1 ] \mathbf{K} = \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} K=fx000fy0cxcy1
    • R: Rectification matrix (stereo cameras only).

      • A rotation matrix aligning the camera coordinate system to the ideal stereo image plane so that epipolar lines in both stereo images are parallel.
      • For monocular cameras R = I \mathbf{R} = \mathbf{I} R=I
    • P: Projection/camera matrix.

      • For monocular cameras
        P = [ f x 0 c x 0 0 f y c y 0 0 0 1 0 ] \mathbf{P} = \begin{bmatrix} f_x & 0 & c_x & 0 \\ 0 & f_y & c_y & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} P=fx000fy0cxcy1000
      • For a stereo pair, the fourth column [Tx Ty 0]’ is related to the position of the optical center of the second camera in the first camera’s frame. We assume Tz = 0 so both cameras are in the same stereo image plane.
        • The first camera
          P = [ f x ′ 0 c x ′ 0 0 f y ′ c y ′ 0 0 0 1 0 ] \mathbf{P} = \begin{bmatrix} f_x' & 0 & c_x' & 0 \\ 0 & f_y' & c_y' & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} P=fx000fy0cxcy1000
        • The second camera
          P = [ f x ′ 0 c x ′ − f x ′ ⋅ B 0 f y ′ c y ′ 0 0 0 1 0 ] \mathbf{P} = \begin{bmatrix} f_x' & 0 & c_x' & -f_x' \cdot B \\ 0 & f_y' & c_y' & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} P=fx000fy0cxcy1fxB00
      • Given a 3D point [ X Y Z ] ′ [X Y Z]' [XYZ], the projection ( x , y ) (x, y) (x,y) of the point onto the rectified image is given by:

      [ u v w ] = P ⋅ [ X Y Z 1 ] , { x = u w y = v w \begin{bmatrix} u \\ v \\ w \end{bmatrix} = \mathbf{P} \cdot \begin{bmatrix} X \\ Y \\ Z \\ 1 \end{bmatrix} , \quad \begin{cases} x = \frac{u}{w} \\ y = \frac{v}{w} \end{cases} uvw=PXYZ1,{x=wuy=wv

    initUndistortRectifyMap

    左右目 分别利用 OpenCV函数 initUndistortRectifyMap 计算 the undistortion and rectification transformation map,得到

    • 左目map: m a p 1 l , m a p 2 l map^l_1, map^l_2 map1l,map2l
    • 右目map: m a p 1 r , m a p 2 r map^r_1, map^r_2 map1r,map2r

    示例代码:

    cv::fisheye::initUndistortRectifyMap(K1, D1, R1, P1, img_size, CV_16SC2, rect_map_[0][0], rect_map_[0][1]);
    cv::fisheye::initUndistortRectifyMap(K2, D2, R2, P2, img_size, CV_16SC2, rect_map_[1][0], rect_map_[1][1]);
    

    Remap

    左右目 分别利用 OpenCV函数 remap 并根据 左右目map 对左右目图像进行 去畸变 和 立体矫正,得到 左右目矫正图像

    示例代码:

    cv::remap(img_l, img_rect_l, rect_map_[0][0], rect_map_[0][1], cv::INTER_LINEAR);
    cv::remap(img_r, img_rect_r, rect_map_[1][0], rect_map_[1][1], cv::INTER_LINEAR);
    

    4. 立体匹配

    根据双目矫正图像,通过 BM或SGM等立体匹配算法 对其进行立体匹配,计算 视差图

    视差计算

    通过 OpenCV函数 stereoBM (block matching algorithm),生成 视差图(Disparity Map) (CV_16S or CV_32F)

    disparity map from stereoBM of OpenCV :
    It has the same size as the input images. When disptype == CV_16S, the map is a 16-bit signed single-channel image, containing disparity values scaled by 16. To get the true disparity values from such fixed-point representation, you will need to divide each disp element by 16. If disptype == CV_32F, the disparity map will already contain the real disparity values on output.

    So if you’ve chosen disptype = CV_16S during computation, you can access a pixel at pixel-position (X,Y) by: short pixVal = disparity.at<short>(Y,X);, while the disparity value is float disparity = pixVal / 16.0f;; if you’ve chosen disptype = CV_32F during computation, you can access the disparity directly: float disparity = disparity.at<float>(Y,X);

    5. 三维重建

    (1)算法1:根据视差图,利用 f ′ f' f B B B 通过几何关系计算 深度值,并利用相机内参计算 三维坐标

    根据上图相似三角形关系,得

    Z B = Z − f B − d w ⟹ Z = B f d w \frac{Z}{B} = \frac{Z-f}{B-d_w} \quad \Longrightarrow \quad Z = \frac{Bf}{d_w} BZ=BdwZfZ=dwBf

    其中, f f f d w d_w dw 分别为 成像平面的焦距和视差,单位均为 物理单位,将其转换为 像素单位,上式写为

    Z = B f ′ d p Z = \frac{B f'}{d_p} Z=dpBf

    其中,

    d p = ( O r − u r ) + ( u l − O l ) = ( u l − u r ) + ( O r − O l ) d_p = (O_r - u_r) + (u_l - O_l) = (u_l - u_r) + (O_r - O_l) dp=(Orur)+(ulOl)=(ulur)+(OrOl)

    最终,深度计算公式如下,通过遍历图像可生成 深度图

    Z = depth = B ⋅ f ′ d p with d p = disp ( u , v ) + ( c x 2 ′ − c x 1 ′ ) Z = \text{depth} = \frac{B \cdot f'}{d_p} \quad \text{with} \quad d_p = \text{disp}(u,v) + ({c_x}_2' - {c_x}_1') Z=depth=dpBfwithdp=disp(u,v)+(cx2cx1)

    根据 小孔成像模型,已知 Z Z Z相机内参 可计算出 三维点坐标,从而可生成 三维点云

    { Z = depth = f ′ ⋅ B d p X = u − c x 1 ′ f ′ ⋅ Z Y = v − c y ′ f ′ ⋅ Z 或 { bd = B d p Z = depth = f ′ ⋅ bd X = ( u − c x 1 ′ ) ⋅ bd Y = ( u − c y ′ ) ⋅ bd \begin{aligned} \begin{cases} Z = \text{depth} = \frac{f' \cdot B}{d_p} \\ X = \frac{u-{c_x}_1'}{f'} \cdot Z \\ Y = \frac{v-{c_y}'}{f'} \cdot Z \end{cases} \end{aligned} \text{或} \begin{aligned} \begin{cases} \text{bd} = \frac{B}{d_p}\\ Z = \text{depth} = f' \cdot \text{bd} \\ X = (u-{c_x}_1') \cdot \text{bd} \\ Y = (u-{c_y}') \cdot \text{bd} \end{cases} \end{aligned} Z=depth=dpfBX=fucx1ZY=fvcyZbd=dpBZ=depth=fbdX=(ucx1)bdY=(ucy)bd

    其中, disp ( u , v ) \text{disp}(u,v) disp(u,v) 代表 视差图 坐标值

    (2)算法2:根据视差图,利用 Q Q Q 矩阵 计算 三维点坐标(reprojectImageTo3D

    [ X ′ Y ′ Z ′ W ] = Q ⋅ [ u v disp ( u , v ) 1 ] \begin{bmatrix} X' \\ Y' \\ Z' \\ W \end{bmatrix} = Q \cdot \begin{bmatrix} u \\ v \\ \text{disp}(u,v) \\ 1 \end{bmatrix} XYZW=Quvdisp(u,v)1

    最终,三维点坐标为

    [ X Y Z ] = [ X ′ W Y ′ W Z ′ W ] \begin{bmatrix} X \\ Y \\ Z \end{bmatrix} = \begin{bmatrix} \frac{X'}{W} \\[2ex] \frac{Y'}{W} \\[2ex] \frac{Z'}{W} \end{bmatrix} XYZ=WXWYWZ

    深度图 图像类型

    • 单位meter --> 32FC1
    • 单位millimeter --> 16UC1

    总结

    展开全文
    u011178262 2018-07-22 17:15:08
  • 30.53MB andrea_pirlo 2018-06-11 21:13:55
  • 1.5MB weixin_38731075 2021-03-17 13:41:30
  • 5星
    1017KB leavemyleave 2021-10-15 00:53:40
  • 双目 MATLAB标定 ,查阅博主的【计算机视觉】摄像机标定 matlab toolbox_calib工具箱(单目标定和双目标定) 1 基本原理   得到了立体标定参数之后,就可以把参数放入xml文件,然后用cvLoad读入OpenCV了。具体...

    双目 MATLAB标定 ,查阅博主的【计算机视觉】摄像机标定 matlab toolbox_calib工具箱(单目标定和双目标定)

    1  基本原理

     

    得到了立体标定参数之后,就可以把参数放入xml文件,然后用cvLoad读入OpenCV了。具体的方法可以参照Learning OpenCV第11章的例子,上面就是用cvSave保存标定结果,然后再用cvLoad把之前的标定结果读入矩阵的

    2.14. xml文件示例

    2.14. xml文件示例

     

     

    这里需要注意的是Matlab标定结果中的om向量,这个向量是旋转矩阵通过Rodrigues变换之后得出的结果,如果要在cvStereoRectify中使用的话,需要首先将这个向量用cvRodrigues转换成旋转矩阵。关于Rodrigues变换,Learning OpenCV的第11章也有说明。

     

    2.15. 旋转矩阵的Rodrigues形式表示

    2.15. 旋转矩阵的Rodrigues形式表示

     

     

     

     

     

     

    3. 立体校准和匹配

    有了标定参数,校准的过程就很简单了。

    我使用的是OpenCV中的cvStereoRectify,得出校准参数之后用cvRemap来校准输入的左右图像。这部分的代码参考的是Learning OpenCV 十二章的例子。

     

    校准之后,就可以立体匹配了。立体匹配OpenCV里面有两种方法,一种是Block Matching,一种是Graph Cut。Block Matching用的是SAD方法,速度比较快,但效果一般。Graph Cut可以参考Kolmogrov03的那篇博士论文,效果不错,但是运行速度实在是慢到不能忍。所以还是选择BM。

     

    以下是我用BM进行立体匹配的参数设置

     

     

    [cpp] view plain copy

    1. BMState = cvCreateStereoBMState(CV_STEREO_BM_BASIC,0);  
    2. assert(BMState != 0);  
    3. BMState->preFilterSize=13;  
    4. BMState->preFilterCap=13;  
    5. BMState->SADWindowSize=19;  
    6. BMState->minDisparity=0;  
    7. BMState->numberOfDisparities=unitDisparity*16;  
    8. BMState->textureThreshold=10;  
    9. BMState->uniquenessRatio=20;  
    10. BMState->speckleWindowSize=13;  

     

     

    其中minDisparity这个参数我设置为0是由于我的两个摄像头是前向平行放置,相同的物体在左图中一定比在右图中偏右,如下图3.1所示。所以没有必要设置回搜的参数。

    如果为了追求更大的双目重合区域而将两个摄像头向内偏转的话,这个参数是需要考虑的。

     

    3.1. 校正后的左右视图

    3.1. 校正后的左右视图

     

     

    另外需要提的参数是uniquenessRatio,实验下来,我感觉这个参数对于最后的匹配结果是有很大的影响。uniquenessRatio主要可以防止误匹配,其主要作用从下面三幅图的disparity效果比对就可以看出。在立体匹配中,我们宁愿区域无法匹配,也不要误匹配。如果有误匹配的话,碰到障碍检测这种应用,就会很麻烦。

     

    3.2. UniquenessRatio为0时的匹配图,可以看到大片的误匹配区域

    3.2. UniquenessRatio为0时的匹配图,可以看到大片的误匹配区域

     

    3.3. UniquenessRatio为10时的disparity map, 可以看到误匹配被大量减少了

    3.3. UniquenessRatio为10时的disparity map, 可以看到误匹配被大量减少了, 但还是有噪点

     

    3.4. UniquenessRatio为20时的disparity map, 可以看到误匹配基本被去除了, 点云干净了很多

    3.4. UniquenessRatio为20时的disparity map, 可以看到误匹配基本被去除了, 点云干净了很多

     

     

     

    关于cvFindStereoCorrespondenceBM这个函数的源代码,曾经做过比较详细的研究,过一段时间也会把之前写的代码注释整理一下,发篇博文。

     

     

     

     

    4. 实际距离的测量

    在用cvFindStereoCorrespondenceBM得出disparity map之后,还需要通过cvReprojectImageTo3D这个函数将单通道Disparity Map转换成三通道的实际坐标矩阵。

    具体的数学原理可以参考下面这个公式(from chenyusiyuan http://blog.csdn.net/chenyusiyuan/archive/2009/12/25/5072597.aspx,实际深度的一些问题这篇博文中也有提到)

    距离转换公式

    4.1 距离转换公式

     

     

    但是在实际操作过程中,用cvReprojectImageTo3D得到的数据并未如实际所想,生成深度矩阵所定义的世界坐标系我就一直没弄清楚。这在下面的例子中会详细说明,希望这方面的专家能帮忙解答一下:

     

     

    图4.2是测量时的实际场景图,场景中主要测量的三个物体就是最前面的利乐包装盒、中间的纸杯、和最远的塑料瓶。

     

     

    4.2. 实际场景中三个待测物体的位置

    4.2. 实际场景中三个待测物体的位置

     

     

    图4.3是校准后的左右图和匹配出来的disparity map,disparity窗口中是实际的点云,object窗口是给disparity map加了个阈值之后得到的二值图,主要是为了分割前景和背景。可以看到要测的三个物体基本被正确地分割出来了

     

    4.3. 双目摄像头得到的disparity map

    4.3. 双目摄像头得到的disparity map

     

    图4.4是在disparity窗口中选取一个点后然后在实际坐标矩阵中得到的对应三维信息,在这里,我在三个物体的点云上各选一个点来代表一个物体实际的坐标信息。(这里通过鼠标获取一点坐标信息的方法参考的是opencv sample里的watershed.cpp)

     

     

    4.4. 对应点的三维坐标

    4.4. 对应点的三维坐标

     

     

    在这里可以看到,(265, 156)也就是利乐包装盒的坐标是(13, 12, -157),(137, 142)纸杯的坐标是(77, 30, -312),(95, 115)塑料瓶的坐标是(144, 63, -482)。

    补充一下:为了方便显示,所以视差图出来之后进行了一个0-255的normalize,所以value值的前一个是normalize之后点的灰度值,后一个是normalize之前点的实际视差图。

    由cvFindStereoCorrespondenceBM算法的源代码:

     

    dptr[y*dstep] = (short)(((ndisp - mind - 1 + mindisp)*256 + (d != 0 ? (p-n)*128/d : 0) + 15) >> 4);

    其中

    ndisp是ndisp = state->numberOfDisparities;

    mindisp是mindisp = state->minDisparity;

    mind就是sad得出的视差

    实际视差大约是(64-mind-1)*256=1163, 基本是对的, 后面一项修正值在通常情况下可以忽略

     

     

     

    目前我还是不是很清楚立体坐标系原点(左相机坐标原点)和尺度,但是从这三个点的z坐标可以大致看出这三个物体的距离差大概是1:2:3,基本与实际场景中物体的位置一致。因此,可以通过这种方法确定出物体的大致距离信息。

     

    但是,如果就从摄像头参数本身来测量距离的话,就不是很明白了,还求这方面的大牛解答。

     

    5.  一些问题

    5.1 关于StereoCalibrate

    OpenCV自带的cvStereoCalibrate感觉不怎么好用,用这个函数求出的内参外参和旋转平移矩阵进行校准,往往无法达到行对准,有时甚至会出现比较可怕的畸变。在看了piao的http://www.opencv.org.cn/forum/viewtopic.php?f=1&t=4603帖子之后,也曾经尝试过现用cvCalibrateCamera2单独标定(左右各20幅图),得出的结果基本和Matlab单独标定的相同,然后再在cvStereoCalibrate中将参数设成CV_CALIB_USE_INTRINSIC_GUESS,用来细化内参数和畸变参数,结果得出的标定结果就又走样了。

    不知道有谁在这方面有过成功经验的,可以出来分享一下。毕竟用Matlab工具箱还是麻烦了些。

     

     

     

    5.2 Translation向量以及立体匹配得出的世界坐标系

    Learning OpenCV中对于Translation和Rotation的图示是这样的

    5.1. Learning OpenCV中的图示

    5.1. Learning OpenCV中的图示

     

    可是在实验过程中发现,如果将Translation向量按尺度缩放,对于StereoRectify之后的左右视图不会有变化,比如将T = [ -226.73817   -0.62302  8.93984 ] ,变成T = [ -22.673817   -0.062302  0.893984 ],在OpenCV中显示的结果不会有任何变化。而且我如果修改其中的一个参量的话,左右视图发生的变化也不是图5.1中所示的那种变化(比如把x缩小,那么视图发生的变化不是往x轴方向的平移)。

     

    因此又回到了老问题,这里这些坐标的尺度究竟是什么?通过ReprojectTo3D那个函数得到的三维坐标又是以哪个点为原点,那三个方向为x,y,z轴的? 

     

    补充: 对这个问题的解答来自于和maxwellsdemon的讨论

    他的解释如下:rotation是两者的旋转角度的关系,但是你要把它矫正平行,也是需要translation matrix的。你可以设想,两个看似已经平行了的摄像头,但是深度上放置的有差距,那么在矫正的时候会议translation matrix所对应的角度或者直线为基准,二者旋转一个小角度,使得完全平行。

     

     

     

     

     

    目前立体匹配算法是计算机视觉中的一个难点和热点,算法很多,但是一般的步骤是:

     

     

    A、匹配代价计算

    匹配代价计算是整个立体匹配算法的基础,实际是对不同视差下进行灰度相似性测量。常见的方法有灰度差的平方SD(squared intensity differences),灰度差的绝对值AD(absolute intensity differences)等。另外,在求原始匹配代价时可以设定一个上限值,来减弱叠加过程中的误匹配的影响。以AD法求匹配代价为例,可用下式进行计算,其中T为设定的阈值。

     

    图18

     

     

    B、 匹配代价叠加

    一般来说,全局算法基于原始匹配代价进行后续算法计算。而区域算法则需要通过窗口叠加来增强匹配代价的可靠性,根据原始匹配代价不同,可分为:

     

    图19

     

     

    C、 视差获取

    对于区域算法来说,在完成匹配代价的叠加以后,视差的获取就很容易了,只需在一定范围内选取叠加匹配代价最优的点(SAD和SSD取最小值,NCC取最大值)作为对应匹配点,如胜者为王算法WTA(Winner-take-all)。而全局算法则直接对原始匹配代价进行处理,一般会先给出一个能量评价函数,然后通过不同的优化算法来求得能量的最小值,同时每个点的视差值也就计算出来了。

     

     

    D、视差细化(亚像素级)

    大多数立体匹配算法计算出来的视差都是一些离散的特定整数值,可满足一般应用的精度要求。但在一些精度要求比较高的场合,如精确的三维重构中,就需要在初始视差获取后采用一些措施对视差进行细化,如匹配代价的曲线拟合、图像滤波、图像分割等。

     

    有关立体匹配的介绍和常见匹配算法的比较,推荐大家看看Stefano Mattoccia 的讲义 Stereo Vision: algorithms and applications,190页的ppt,讲解得非常形象详尽。

     

     

     

     

     

    1. opencv2.1和opencv2.0在做stereo vision方面有什么区别了?

    2.1版增强了Stereo Vision方面的功能:

    (1) 新增了 SGBM 立体匹配算法(源自Heiko Hirschmuller的《Stereo Processing by Semi-global Matching and Mutual Information》),可以获得比 BM 算法物体轮廓更清晰的视差图(但低纹理区域容易出现横/斜纹路,在 GCstate->fullDP 选项使能时可消减这种异常纹路,但对应区域视差变为0,且运行速度会有所下降),速度比 BM 稍慢, 352*288的帧处理速度大约是 5 帧/秒;

    (2) 视差效果:BM < SGBM < GC;处理速度:BM > SGBM > GC ;

    (3) BM 算法比2.0版性能有所提升,其状态参数新增了对左右视图感兴趣区域 ROI 的支持(roi1 和 roi2,由stereoRectify函数产生);

    (4) BM 算法和 GC 算法的核心代码改动不大,主要是面向多线程运算方面的(由 OpenMP 转向 Intel TBB);

    (5) cvFindStereoCorrespondenceBM 函数的disparity参数的数据格式新增了 CV_32F 的支持,这种格式的数据给出实际视差,而 2.0 版只支持 CV_16S,需要除以 16.0 才能得到实际的视差数值。

     

     

    2. 用于立体匹配的图像可以是彩色的吗?

    在OpenCV2.1中,BM和GC算法只能对8位灰度图像计算视差,SGBM算法则可以处理24位(8bits*3)彩色图像。所以在读入图像时,应该根据采用的算法来处理图像:

     

     

    int color_mode = alg == STEREO_SGBM ? 1 : 0;
    //
    // 载入图像
    cvGrabFrame( lfCam );
    cvGrabFrame( riCam );
    frame1 = cvRetrieveFrame( lfCam );
    frame2 = cvRetrieveFrame( riCam );
    if(frame1.empty()) break;
    resize(frame1, img1, img_size, 0, 0);
    resize(frame2, img2, img_size, 0, 0);
    // 选择彩色或灰度格式作为双目匹配的处理图像
    if (!color_mode && cn>1)
    {
    cvtColor(img1, img1gray, CV_BGR2GRAY);
    cvtColor(img2, img2gray, CV_BGR2GRAY);
    img1p = img1gray;
    img2p = img2gray;
    }
    else
    {
    img1p = img1;
    img2p = img2;
    }

     

     

    3. 怎样获取与原图像有效像素区域相同的视差图?

    OpenCV2.0及以前的版本中,所获取的视差图总是在左侧和右侧有明显的黑色区域,这些区域没有有效的视差数据。视差图有效像素区域与视差窗口(ndisp,一般取正值且能被16整除)和最小视差值(mindisp,一般取0或负值)相关,视差窗口越大,视差图左侧的黑色区域越大,最小视差值越小,视差图右侧的黑色区域越大。其原因是为了保证参考图像(一般是左视图)的像素点能在目标图像(右视图)中按照设定的视差匹配窗口匹配对应点,OpenCV 只从参考图像的第 (ndisp - 1 + mindisp) 列开始向右计算视差,第 0 列到第 (ndisp - 1 + mindisp) 列的区域视差统一设置为 (mindisp - 1) *16;视差计算到第 width + mindisp 列时停止,余下的右侧区域视差值也统一设置为 (mindisp - 1) *16。  

    00177 static const int DISPARITY_SHIFT = 4;
    …
    00411     int ndisp = state->numberOfDisparities;
    00412     int mindisp = state->minDisparity;
    00413     int lofs = MAX(ndisp - 1 + mindisp, 0);
    00414     int rofs = -MIN(ndisp - 1 + mindisp, 0);
    00415     int width = left->cols, height = left->rows;
    00416     int width1 = width - rofs - ndisp + 1;
    …
    00420     short FILTERED = (short)((mindisp - 1) << DISPARITY_SHIFT);
    …
    00466     // initialize the left and right borders of the disparity map
    00467     for( y = 0; y < height; y++ )
    00468     {
    00469         for( x = 0; x < lofs; x++ )
    00470             dptr[y*dstep + x] = FILTERED;
    00471         for( x = lofs + width1; x < width; x++ )
    00472             dptr[y*dstep + x] = FILTERED;
    00473     }
    00474     dptr += lofs;
    00475
    00476     for( x = 0; x < width1; x++, dptr++ )
    
    …
    

     

    这样的设置很明显是不符合实际应用的需求的,它相当于把摄像头的视场范围缩窄了。因此,OpenCV2.1 做了明显的改进,不再要求左右视图和视差图的大小(size)一致,允许对视差图进行左右边界延拓,这样,虽然计算视差时还是按上面的代码思路来处理左右边界,但是视差图的边界得到延拓后,有效视差的范围就能够与对应视图完全对应。具体的实现代码范例如下:

     

    //
    // 对左右视图的左边进行边界延拓,以获取与原始视图相同大小的有效视差区域
    copyMakeBorder(img1r, img1b, 0, 0, m_nMaxDisp, 0, IPL_BORDER_REPLICATE);
    copyMakeBorder(img2r, img2b, 0, 0, m_nMaxDisp, 0, IPL_BORDER_REPLICATE);
    
    //
    // 计算视差
    if( alg == STEREO_BM )
    {
    	bm(img1b, img2b, dispb);
    	// 截取与原始画面对应的视差区域(舍去加宽的部分)
    	displf = dispb.colRange(m_nMaxDisp, img1b.cols);	
    }
    else if(alg == STEREO_SGBM)
    {
    	sgbm(img1b, img2b, dispb);
    	displf = dispb.colRange(m_nMaxDisp, img1b.cols);
    }
    

     

     

     

     

    4. cvFindStereoCorrespondenceBM的输出结果好像不是以像素点为单位的视差?

    @scyscyao:在OpenCV2.0中,BM函数得出的结果是以16位符号数的形式的存储的,出于精度需要,所有的视差在输出时都扩大了16倍(2^4)。其具体代码表示如下:

    dptr[y*dstep] = (short)(((ndisp - mind - 1 + mindisp)*256 + (d != 0 ? (p-n)*128/d : 0) + 15) >> 4);

    可以看到,原始视差在左移8位(256)并且加上一个修正值之后又右移了4位,最终的结果就是左移4位。

    因此,在实际求距离时,cvReprojectTo3D出来的X/W,Y/W,Z/W都要乘以16 (也就是W除以16),才能得到正确的三维坐标信息。”

     

    OpenCV2.1中,BM算法可以用 CV_16S 或者 CV_32F 的方式输出视差数据,使用32位float格式可以得到真实的视差值,而CV_16S 格式得到的视差矩阵则需要 除以16 才能得到正确的视差。另外,OpenCV2.1另外两种立体匹配算法 SGBM 和 GC 只支持 CV_16S 格式的 disparity 矩阵

     

     

     

    5. 如何设置BM、SGBM和GC算法的状态参数?

     

    (1)StereoBMState

    // 预处理滤波参数

    • preFilterType:预处理滤波器的类型,主要是用于降低亮度失真(photometric distortions)、消除噪声和增强纹理等, 有两种可选类型:CV_STEREO_BM_NORMALIZED_RESPONSE(归一化响应) 或者 CV_STEREO_BM_XSOBEL(水平方向Sobel算子,默认类型), 该参数为 int 型;
    • preFilterSize:预处理滤波器窗口大小,容许范围是[5,255],一般应该在 5x5..21x21 之间,参数必须为奇数值, int 型
    • preFilterCap:预处理滤波器的截断值,预处理的输出值仅保留[-preFilterCap, preFilterCap]范围内的值,参数范围:1 - 31(文档中是31,但代码中是 63), int

    // SAD 参数

    • SADWindowSize:SAD窗口大小,容许范围是[5,255],一般应该在 5x5 至 21x21 之间,参数必须是奇数,int 型
    • minDisparity:最小视差默认值为 0, 可以是负值,int 型
    • numberOfDisparities:视差窗口,即最大视差值与最小视差值之差, 窗口大小必须是 16 的整数倍,int 型

    // 后处理参数

    • textureThreshold:低纹理区域的判断阈值。如果当前SAD窗口内所有邻居像素点的x导数绝对值之和小于指定阈值,则该窗口对应的像素点的视差值为 0(That is, if the sum of absolute values of x-derivatives computed over SADWindowSize by SADWindowSize pixel neighborhood is smaller than the parameter, no disparity is computed at the pixel),该参数不能为负值,int 型
    • uniquenessRatio:视差唯一性百分比, 视差窗口范围内最低代价是次低代价的(1 + uniquenessRatio/100)倍时,最低代价对应的视差值才是该像素点的视差,否则该像素点的视差为 0 (the minimum margin in percents between the best (minimum) cost function value and the second best value to accept the computed disparity, that is, accept the computed disparity d^ only if SAD(d) >= SAD(d^) x (1 + uniquenessRatio/100.) for any d != d*+/-1 within the search range ),该参数不能为负值,一般5-15左右的值比较合适,int 型
    • speckleWindowSize:检查视差连通区域变化度的窗口大小, 值为 0 时取消 speckle 检查,int 型
    • speckleRange:视差变化阈值,当窗口内视差变化大于阈值时,该窗口内的视差清零,int 型

    // OpenCV2.1 新增的状态参数

    • roi1, roi2:左右视图的有效像素区域,一般由双目校正阶段的 cvStereoRectify 函数传递,也可以自行设定。一旦在状态参数中设定了 roi1 和 roi2,OpenCV 会通过cvGetValidDisparityROI 函数计算出视差图的有效区域,在有效区域外的视差值将被清零。
    • disp12MaxDiff:左视差图(直接计算得出)和右视差图(通过cvValidateDisparity计算得出)之间的最大容许差异。超过该阈值的视差值将被清零。该参数默认为 -1,即不执行左右视差检查。int 型。注意在程序调试阶段最好保持该值为 -1,以便查看不同视差窗口生成的视差效果。具体请参见《使用OpenGL动态显示双目视觉三维重构效果示例》一文中的讨论。

    在上述参数中,对视差生成效果影响较大的主要参数是 SADWindowSize、numberOfDisparities 和 uniquenessRatio 三个,一般只需对这三个参数进行调整,其余参数按默认设置即可

    在OpenCV2.1中,BM算法有C和C++ 两种实现模块。

     

    (2)StereoSGBMState

    SGBM算法的状态参数大部分与BM算法的一致,下面只解释不同的部分:

    • SADWindowSize:SAD窗口大小,容许范围是[1,11],一般应该在 3x3 至 11x11 之间,参数必须是奇数,int 型
    • P1, P2:控制视差变化平滑性的参数。P1、P2的值越大,视差越平滑。P1是相邻像素点视差增/减 1 时的惩罚系数;P2是相邻像素点视差变化值大于1时的惩罚系数。P2必须大于P1。OpenCV2.1提供的例程 stereo_match.cpp 给出了 P1 和 P2 比较合适的数值
    • fullDP:布尔值,当设置为 TRUE 时,运行双通道动态编程算法(full-scale 2-pass dynamic programming algorithm),会占用O(W*H*numDisparities)个字节,对于高分辨率图像将占用较大的内存空间。一般设置为 FALSE

    注意OpenCV2.1的SGBM算法是用C++ 语言编写的,没有C实现模块。与H. Hirschmuller提出的原算法相比,主要有如下变化:

    1. 算法默认运行单通道DP算法,只用了5个方向,而fullDP使能时则使用8个方向(可能需要占用大量内存)。
    2. 算法在计算匹配代价函数时,采用块匹配方法而非像素匹配(不过SADWindowSize=1时就等于像素匹配了)。
    3. 匹配代价的计算采用BT算法("Depth Discontinuities by Pixel-to-Pixel Stereo" by S. Birchfield and C. Tomasi),并没有实现基于互熵信息的匹配代价计算。
    4. 增加了一些BM算法中的预处理和后处理程序。

     

    (3)StereoGCState

    GC算法的状态参数只有两个:numberOfDisparities 和 maxIters ,并且只能通过 cvCreateStereoGCState 在创建算法状态结构体时一次性确定,不能在循环中更新状态信息。GC算法并不是一种实时算法,但可以得到物体轮廓清晰准确的视差图,适用于静态环境物体的深度重构。

    注意GC算法只能在C语言模式下运行,并且不能对视差图进行预先的边界延拓,左右视图和左右视差矩阵的大小必须一致。

     

     

    6. 如何实现视差图的伪彩色显示?

    首先要将16位符号整形的视差矩阵转换为8位无符号整形矩阵,然后按照一定的变换关系进行伪彩色处理。我的实现代码如下: 

     

    // 转换为 CV_8U 格式,彩色显示
    dispLfcv = displf, dispRicv = dispri, disp8cv = disp8;
    if (alg == STEREO_GC)
    {
    	cvNormalize( &dispLfcv, &disp8cv, 0, 256, CV_MINMAX );
    } 
    else
    {
    	displf.convertTo(disp8, CV_8U, 255/(m_nMaxDisp*16.));
    }
    F_Gray2Color(&disp8cv, vdispRGB);
    

     

     

    灰度图转伪彩色图的代码,主要功能是使灰度图中 亮度越高的像素点,在伪彩色图中对应的点越趋向于 红色;亮度越低,则对应的伪彩色越趋向于 蓝色;总体上按照灰度值高低,由红渐变至蓝,中间色为绿色。其对应关系如下图所示:

     

    图20

     

     

    void F_Gray2Color(CvMat* gray_mat, CvMat* color_mat)
    {
    	if(color_mat)
    		cvZero(color_mat);
    		
    	int stype = CV_MAT_TYPE(gray_mat->type), dtype = CV_MAT_TYPE(color_mat->type);
    	int rows = gray_mat->rows, cols = gray_mat->cols;
    
    	// 判断输入的灰度图和输出的伪彩色图是否大小相同、格式是否符合要求
    	if (CV_ARE_SIZES_EQ(gray_mat, color_mat) && stype == CV_8UC1 && dtype == CV_8UC3)
    	{
    		CvMat* red = cvCreateMat(gray_mat->rows, gray_mat->cols, CV_8U);
    		CvMat* green = cvCreateMat(gray_mat->rows, gray_mat->cols, CV_8U);
    		CvMat* blue = cvCreateMat(gray_mat->rows, gray_mat->cols, CV_8U);
    		CvMat* mask = cvCreateMat(gray_mat->rows, gray_mat->cols, CV_8U);
    
    		// 计算各彩色通道的像素值
    		cvSubRS(gray_mat, cvScalar(255), blue);	// blue(I) = 255 - gray(I)
    		cvCopy(gray_mat, red);			// red(I) = gray(I)
    		cvCopy(gray_mat, green);			// green(I) = gray(I),if gray(I) < 128
    		cvCmpS(green, 128, mask, CV_CMP_GE );	// green(I) = 255 - gray(I), if gray(I) >= 128
    		cvSubRS(green, cvScalar(255), green, mask);
    		cvConvertScale(green, green, 2.0, 0.0);
    
    		// 合成伪彩色图
    		cvMerge(blue, green, red, NULL, color_mat);
    
    		cvReleaseMat( &red );
    		cvReleaseMat( &green );
    		cvReleaseMat( &blue );
    		cvReleaseMat( &mask );
    	}
    }
    

     

     

     

    7. 如何将视差数据保存为 txt 数据文件以便在 Matlab 中读取分析?

    由于OpenCV本身只支持 xml、yml 的数据文件读写功能,并且其xml文件与构建网页数据所用的xml文件格式不一致,在Matlab中无法读取。我们可以通过以下方式将视差数据保存为txt文件,再导入到Matlab中。 

     

    void saveDisp(const char* filename, const Mat& mat)		
    {
    	FILE* fp = fopen(filename, "wt");
    	fprintf(fp, "%02d/n", mat.rows);
    	fprintf(fp, "%02d/n", mat.cols);
    	for(int y = 0; y < mat.rows; y++)
    	{
    		for(int x = 0; x < mat.cols; x++)
    		{
    			short disp = mat.at<short>(y, x); // 这里视差矩阵是CV_16S 格式的,故用 short 类型读取
    			fprintf(fp, "%d/n", disp); // 若视差矩阵是 CV_32F 格式,则用 float 类型读取
    		}
    	}
    	fclose(fp);
    }
    

     

    相应的Matlab代码为:

      

    function img = txt2img(filename)
    data = importdata(filename);
    r = data(1);    % 行数
    c = data(2);    % 列数
    disp = data(3:end); % 视差
    vmin = min(disp);
    vmax = max(disp);
    disp = reshape(disp, [c,r])'; % 将列向量形式的 disp 重构为 矩阵形式
    %  OpenCV 是行扫描存储图像,Matlab 是列扫描存储图像
    %  故对 disp 的重新排列是首先变成 c 行 r 列的矩阵,然后再转置回 r 行 c 列
    img = uint8( 255 * ( disp - vmin ) / ( vmax - vmin ) );
    mesh(disp);
    set(gca,'YDir','reverse');  % 通过 mesh 方式绘图时,需倒置 Y 轴方向
    axis tight; % 使坐标轴显示范围与数据范围相贴合,去除空白显示区
    

     

    显示效果如下:

     

    图21

     

     

     

    展开全文
    KYJL888 2017-05-07 21:41:02
  • 634KB weixin_38743481 2019-09-12 15:32:18
  • 基础知识 极线搜索 光度一致性约束 可视性约束 多视角立体技术 基于体素的方法 基于空间patch的方法 深度图融合

    基础知识

      稠密点云重建的目的是在相机位姿已知的前提下,逐像素的计算图像中每一个像素点对应的三维点,得到场景物体表面密集的三维点云。
      稠密点云重建的前提相机的姿态已知,且待重建的三维点具有图像一致性,即三维点位于物体的表面上而非物体内部,那么这个三维点根据相机的内外参投影到不同视角的图像上,以投影点为中心取出两个方形的小patch,那么它们所包含的场景应该是比较类似的.

    极线搜索

      两个视角的图像中,A视角中的一点在不同的深度的三维点投影到B视角中形成一条线,称为极线。当三维点的坐标未知,A图像中的一点,在B图像中可能匹配上对应的极线上的某点,在极线上找到B点的具体位置与相机的姿态之后,可以利用三角量量测的得到三位点的深度信息。

    在这里插入图片描述

      如何判断极线上的某点是否为A点的同名点,主要利用以下约束:光度一致性约束、可视性约束。


    光度一致性约束

      同一个空间的点在不同的视角的投影应当具有相同的光度,重建的核心在于恢复空间中具有光度一致性的点。(要求物体的表面为朗伯面,玻璃、镜子、水面不属于朗伯面)

    光度一致性假设的度量方法
    Sum of squared Differences(SSD):
    ρ S S D ( f , g ) = ∣ ∣ f − g ∣ ∣ 2 \rho_{SSD}(f, g)=||f-g||^{2} ρSSD(f,g)=fg2
    Sum of Absolute Differences(SAD):
    ρ S A D ( f , g ) = ∣ ∣ f − g ∣ ∣ 1 \rho_{SAD}(f,g)=||f-g||_{1} ρSAD(f,g)=fg1
    Normalized Cross Correlation(NCC):
    ρ N C C ( f , g ) = ( f − f ‾ ) ( g − g ‾ ) δ f δ g \rho_{NCC}(f,g)=\frac{(f-\overline f)(g-\overline g)}{\delta_{f}\delta_{g}} ρNCC(f,g)=δfδg(ff)(gg)
    其中f代表以A视图中某点的领域的像素组成的向量,g代表B视图中对应同名点邻域像素组成的向量


    可视性约束

      1、图像上的点不能被遮挡。 2、重建的点前面不能出现点。 3、不能出现在物体的内部。
    在这里插入图片描述


    多视角立体技术


    基于体素的方法

    规则划分

      基于体素的方法等价于3D空间中的像素标定问题,在物体内部的点标记为1,在物体外部的点,标记为0,介于0-1之间的点就是物体的表面。
      标记完之后,使用光度一致性约束和可视性约束来对物体表面的点进行计算,评估标定质量。基于体素的方法的缺点是计算量大,在分辨率很大的时候,计算量会增长很多。

    在这里插入图片描述

    不规则划分

    在点云较为稀疏的区域采用分辨率较小的体素,在点云较为稠密与细节比较精确的的区域采用分辨率较高的体素,可以保证物体重建的精度,另一方面可以自适应调整四面体的大小,减少计算量。

    在这里插入图片描述

    基于体素的方法的优化

    典型的MRF离散化优化问题

    参考图像上的每个像素都分配一个标签(内部或者外部)(无向图的优化)
    E ( f ) = E d a t a ( f ) + E s m o o t h n e s s ( f ) + E v i s i b i l i t y ( f ) E(f)=E_{data}(f) + E_{smoothness}(f) + E_{visibility}(f) E(f)=Edata(f)+Esmoothness(f)+Evisibility(f)
    其中 E ( f ) E(f) E(f)代表能量模型。
    f表示无向图的所有顶点的向量,需要对每个顶点赋值
    E d a t a E_{data} Edata是光度一致性假设,在不同视图上观察到该点的纹理或者光度尽量相似。
    E s m o o t h n e s s E_{smoothness} Esmoothness平滑项约束,相邻的体素(四面体)的标签趋于一致。
    E v i s i b i l i t y E_{visibility} Evisibility不能被别的东西遮挡,不能遮挡别的体素。

    基于体素的方法的优小结

    在这里插入图片描述

    优点:1、生成规则的点云
       2、便于提取物体的平面

    缺点:1、精度受到空间划分分辨率的影响
       2、难以处理精度高、规模大的场景


    基于空间patch扩散的方法(pmvs)

    1、假设空间中的3d矩形patch
    2、通过一定规则的扩张方法,使得pathch覆盖物体表面

    在这里插入图片描述

    过程:

    1. 初始种子点生成采用sift、HOG等稀疏特征点,对这些点进行深度和法向量的重建
    2. 扩张过程对已重建三维点的邻域进行匹配
    3. 滤波过程采用光度一致性约束和可视性约束来去除噪点

    在这里插入图片描述

    3D patch定义:
    位置、法向量(每个patch具有25个点, pathch的中心为点的位置,法向量为三维点邻域的法向量)

    在这里插入图片描述
    在这里插入图片描述

    对三维空间上的一个pathch投影到不同视角中,计算NCC
    1.在图像上均匀计算HOG/Harris特征
    2.沿极线进行搜索找到匹配特征点
    3.对匹配对进行三角化建立patch,
      中心: 三角化确定
      法向量: 指向参考图像
      可视图像 V ( p ) V(p) V(p): 法向量和视角的夹角足够小,通过光度一致性约束对可视图像进行筛选 V ∗ ( p ) = { I ∣ I ∈ V ( P ) , h ( p , I , R ( P ) ) ≤ a } V^{*}(p)=\{I|I\in V(P), h(p,I,R(P))\le a\} V(p)={IIV(P),h(p,I,R(P))a}
    4.对patch位置和法向量进行优化 m a x c ( p ) , n ( p ) g ∗ ( p ) = 1 ∣ V ∗ ( p ) \ R ( p ) ∣ ∑ I ∈ V ∗ ( p ) \ R ( p ) h ( p , I , R ( p ) ) \mathop{max}\limits_{c(p), n(p)} g^{*}(p) = \frac{1}{|V^{*}(p) \backslash R(p)|}\sum_{I\in V^{*}(p) \backslash R(p)}h(p,I,R(p)) c(p),n(p)maxg(p)=V(p)\R(p)1IV(p)\R(p)h(p,I,R(p))

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述


    深度图融合(重点)

    基本步骤:

    1. 为图像选择邻域图像构成立体图像对(通过图像pair进行深度图重建)
    2. 计算每一幅参考图像的深度图,得到每幅参考图像的点云
    3. 对所有深度图进行融合

    全局视角的选择(建立候选视角池)

    初步筛选出一些适合用于重建的图像,局部视角从全局视角中选择.
    1、图像具有相同的内容、外观和尺度
    2、图像具有较大的视差(宽基线)

    在这里插入图片描述
    在这里插入图片描述

    其中, F V F_{V} FV表示V视角, F R F_{R} FR表示R视角, w N ( f ) w_{N}(f) wN(f)代表相机的三角测量夹角大小, w s ( f ) w_{s}(f) ws(f)代表两幅图像的尺度相似性。

    图像分辨率的估计

    用于衡量图像的分辨率,图像的一个像素的宽度对应的三维空间中的物体尺寸。

    在这里插入图片描述

    相关视角的选择(针对每个像素)

      全局视角的选择时固定的,局部视角的选择是针对每个像素的,每个像素在重建过程中不需要用到全部的全局视角,有些视角存在遮挡。局部视角选择是从全局视角中再一次进行筛选得到。一般从全局视角中选择四个当作局部视角进行优化。

    局部视角的选择:
    1、NCC值确定候选视角
    2、实现要足够分散(不共面)
    l R ( V ) = g R ∗ ∏ V ′ ∈ A w e ( V , V ′ ) l_{R}(V)=g_{R}*\prod_{V^{'}\in A} w_{e}(V, V^{'}) lR(V)=gRVAwe(V,V)

    深度图融合算法整体思路

    区域生长法扩张(宽度优先搜索):

    1、根据重建的置信度(NCC)建立优先级队列(排序)
    2、从初始的稀疏特征点开始深度估计
    3、对每个种子点进行非线性深度优化(核心
    4、每次优化完成判断两种情况,满足的话将邻域像素添加到队列中当作种子点,并将当前像素patch的信息作为邻域像素的初始值:(1)邻域没有深度值(2)当前像素的置信度高于邻域像素一定范围(当前点因为重建优化完,置信度得到提升,和邻域拉开很大差距,说明这个区域很需要优化,有很大的优化空间

    在这里插入图片描述

      新建patch的可视图像和法向量的初始值等同于邻域patch
    n ( p ′ ) = n ( p ) V ( p ′ ) = V ( P ) n(p^{'})=n(p)\\V(p^{'})=V(P) n(p)=n(p)V(p)=V(P)

    深度值非线性优化

      假设参考视角为 R , k ∈ N R, k\in N R,kN为与参考视角相关的局部邻域视角, P k P_{k} Pk为三维点投影到第k个视角的投影矩阵. 每个像素提供了一个初始的深度和法向量和局部视角后,对深度值进行非线性优化。建立几何模型如下:

    在这里插入图片描述

      其中, o R o_{R} oR表示参考视角的相机中心,(s, t)表示参考视角像素点,其初始深度为 h ( s , t ) h(s,t) h(s,t),已知相机姿态, 则对应的三维空间中的射线的单位向量为 r → R ( s , t ) \overrightarrow{r}_{R}(s, t) r R(s,t), 则该像素对应的三维空间的点为:
    X R ( s , t ) = O R + h ( s , t ) r → R ( s , t ) X_{R}(s, t)=O_{R}+h(s,t)\overrightarrow{r}_{R}(s, t) XR(s,t)=OR+h(s,t)r R(s,t)

    为了对像素点 ( s , t ) (s,t) (s,t)的深度和三维中的法向量进行优化,以该像素为中心,建立nxn的patch, pathch对应三维中间中的一个很小的平面.同时引入两个变量来帮助表示patch中每个像素的三维坐标. 则像素 ( s + i , t + j ) (s+i, t+j) (s+i,t+j)对应的深度为:
    h ( s + i , t + j ) = h ( s , t ) + i ∗ h s ( s , t ) + j ∗ h t ( s , t ) h(s+i, t+j)=h(s,t)+i*h_{s}(s,t)+j*h_{t}(s,t) h(s+i,t+j)=h(s,t)+ihs(s,t)+jht(s,t)

    假设 ( s + i , t + j ) (s+i, t+j) (s+i,t+j)处的射线方向近似为 r → R ( s , t ) \overrightarrow{r}_{R}(s, t) r R(s,t),则像素 ( s + i , t + j ) (s+i, t+j) (s+i,t+j)对应的三维坐标为:
    X R ( i , j ) = O R + r → R ( s , t ) ( h ( s , t ) + i ∗ h s + j ∗ h t ) X_{R}(i,j)=O_{R}+\overrightarrow{r}_{R}(s, t)(h(s,t)+i*h_{s}+j*h_{t}) XR(i,j)=OR+r R(s,t)(h(s,t)+ihs+jht)

    光度模型

      对朗伯反射进行简单假设,以提升颜色匹配的准确度. 为邻域内的每个视角K分配一个颜色尺度 c k = [ c k r , c k g , c k b ] T ∈ R 3 × 1 c_{k}=[c_{k}^{r}, c_{k}^{g}, c_{k}^{b}]^{T} \in{R^{3\times 1}} ck=[ckr,ckg,ckb]TR3×1,如果深度 h ( s , t ) h(s,t) h(s,t)估计准确,则有:
    I R ( s + i , t + j ) = c k ⋅ I k ( P k ( X R ( s + i , t + j ) ) ) I_{R}(s+i,t+j)=c_{k}\centerdot I_{k}(P_{k}(X_{R}(s+i, t+j))) IR(s+i,t+j)=ckIk(Pk(XR(s+i,t+j)))

    对所有的邻域视角和所有的patch中的patch中的三维点成立,其中 I R ( . ) I_{R}(.) IR(.)表示取参考视角对应像素处的颜色(3通道,为向量形式), I k ( . ) I_{k}(.) Ik(.)表示取第k个视角中图像上对应像素处的颜色,表示参考视角的一点和其他视角在对应同一点上的颜色一致.

    结合几何模型和光度模型,可以得到深度值和法向量的优化数学模型如下:
    E = ∑ k ∈ N ∑ i = − n − 1 2 n − 1 2 ∑ j = − n − 1 2 n − 1 2 [ I R ( s + i , t + j ) − c k . I k ( P k ( X R ( s + i , t + j ) ) ) ] 2 E=\sum_{k\in N} \sum_{i=-\frac{n-1}{2}}^{\frac{n-1}{2}} \sum_{j=-\frac{n-1}{2}}^{\frac{n-1}{2}}[I_{R}(s+i,t+j)-c_{k}.I_{k}(P_{k}(X_{R}(s+i, t+j)))]^{2} E=kNi=2n12n1j=2n12n1[IR(s+i,t+j)ck.Ik(Pk(XR(s+i,t+j)))]2
      其中 X R ( s + i , t + j ) X_{R}(s+i, t+j) XR(s+i,t+j)表示将参考图像中的 ( s + i , t + k ) (s+i, t+k) (s+i,t+k)映射到三维空间中的坐标,乘以 P k P_{k} Pk表示将该点的三维坐标投影到第k个视角图像上的像素坐标,对这个像素坐标取 I k I_{k} Ik表示将第k个视角图像上这个点的像素颜色值取出来,再乘以前面的 c k c_{k} ck,表示一个调节因子,增加鲁棒性

      为了书写方便 I R ( s + i , t + j ) I_{R}(s+i,t+j) IR(s+i,t+j)简记为 I R ( i , j ) I_{R}(i,j) IR(i,j), 将 I k ( P k ( X ( s + i , t + j ) ) ) I_{k}(P_{k}(X(s+i,t+j))) Ik(Pk(X(s+i,t+j)))简记为 I k ( i , j ) I_{k}(i,j) Ik(i,j),将 h s ( s , t ) , h t ( s , t ) h_{s}(s,t), h_{t}(s,t) hs(s,t),ht(s,t)分别简记为 h s , h t h_{s}, h_{t} hs,ht于是,上式可以简记为:
    E = ∑ i j k [ I R ( i , j ) − c k . I k ( i , j ) ] 2 E=\sum_{ijk}[I_{R}(i,j)-c_{k}.I_{k}(i,j)]^{2} E=ijk[IR(i,j)ck.Ik(i,j)]2
    i , j ∈ [ − n − 1 2 , n − 1 2 ] i,j\in[-\frac{n-1}{2}, \frac{n-1}{2}] i,j[2n1,2n1]是patch中的采样点,k表示视角的个数,需要优化的变量为三维坐标的深度 h ( s , t ) , h s , h t h(s,t), h_{s},h_{t} h(s,t),hs,ht和颜色尺度因子 { c k } \{c_{k}\} {ck}

    颜色尺度的优化:

    ∂ E ∂ c k r = ∑ i j r k r ( i , j ) I R r ( i , j ) + c k r ∑ i j ( I k r ( i , j ) ) 2 ∂ E ∂ c k g = ∑ i j r k g ( i , j ) I R g ( i , j ) + c k g ∑ i j ( I k g ( i , j ) ) 2 ∂ E ∂ c k b = ∑ i j r k b ( i , j ) I R r ( i , j ) + c k b ∑ i j ( I k b ( i , j ) ) 2 c k = [ ∑ i j I K r ( i , j ) I R r ( i , j ) ∑ i j ( I K r ( i , j ) ) 2 , ∑ i j I K g ( i , j ) I R g ( i , j ) ∑ i j ( I K g ( i , j ) ) 2 , ∑ i j I K b ( i , j ) I R b ( i , j ) ∑ i j ( I K b ( i , j ) ) 2 ] T \frac{\partial E}{\partial c_{k}^{r}} = \sum_{ij}r_{k}^{r}(i,j)I_{R}^{r}(i,j) + c_{k}^{r}\sum_{ij}(I_{k}^{r}(i,j))^{2} \\\frac{\partial E}{\partial c_{k}^{g}} = \sum_{ij}r_{k}^{g}(i,j)I_{R}^{g}(i,j) + c_{k}^{g}\sum_{ij}(I_{k}^{g}(i,j))^{2} \\\frac{\partial E}{\partial c_{k}^{b}} = \sum_{ij}r_{k}^{b}(i,j)I_{R}^{r}(i,j) + c_{k}^{b}\sum_{ij}(I_{k}^{b}(i,j))^{2} \\c_{k}=[\frac{\sum_{ij}I_{K}^{r}(i,j)I_{R}^{r}(i,j)}{\sum_{ij}(I_{K}^{r}(i,j))^2}, \frac{\sum_{ij}I_{K}^{g}(i,j)I_{R}^{g}(i,j)}{\sum_{ij}(I_{K}^{g}(i,j))^2}, \frac{\sum_{ij}I_{K}^{b}(i,j)I_{R}^{b}(i,j)}{\sum_{ij}(I_{K}^{b}(i,j))^2}]^{T} ckrE=ijrkr(i,j)IRr(i,j)+ckrij(Ikr(i,j))2ckgE=ijrkg(i,j)IRg(i,j)+ckgij(Ikg(i,j))2ckbE=ijrkb(i,j)IRr(i,j)+ckbij(Ikb(i,j))2ck=[ij(IKr(i,j))2ijIKr(i,j)IRr(i,j),ij(IKg(i,j))2ijIKg(i,j)IRg(i,j),ij(IKb(i,j))2ijIKb(i,j)IRb(i,j)]T

    颜色尺度的优化:

      分别对上式求三个偏导,并令结果为0,得:
    c k r = ∑ i j I K r ( i , j ) I R r ( i , j ) ∑ i j ( I K r ( i , j ) ) 2 c k g = ∑ i j I K g ( i , j ) I R g ( i , j ) ∑ i j ( I K g ( i , j ) ) 2 c k b = ∑ i j I K b ( i , j ) I R b ( i , j ) ∑ i j ( I K b ( i , j ) ) 2 ] T c_{k}^{r} = \frac{\sum_{ij}I_{K}^{r}(i,j)I_{R}^{r}(i,j)}{\sum_{ij}(I_{K}^{r}(i,j))^2} \\c_{k}^{g}= \frac{\sum_{ij}I_{K}^{g}(i,j)I_{R}^{g}(i,j)}{\sum_{ij}(I_{K}^{g}(i,j))^2} \\c_{k}^{b}=\frac{\sum_{ij}I_{K}^{b}(i,j)I_{R}^{b}(i,j)}{\sum_{ij}(I_{K}^{b}(i,j))^2}]^{T} ckr=ij(IKr(i,j))2ijIKr(i,j)IRr(i,j)ckg=ij(IKg(i,j))2ijIKg(i,j)IRg(i,j)ckb=ij(IKb(i,j))2ijIKb(i,j)IRb(i,j)]T

    h ( s , t ) , h s , h t h(s,t), h_{s},h_{t} h(s,t),hs,ht的优化:

      待优化的模型为最小二乘的非线性优化问题,无法得到该函数的解析解,采用迭代的方式求取其近似解。首先引入 I R ( i , j ) I_{R}(i,j) IR(i,j)关于 h ( s , t ) , h s , h t h(s, t), h_{s}, h_{t} h(s,t),hs,ht的线性表达式:
    I R ( i , j ) ≈ c k ⋅ I k ( P k ( O R + r → R ( s , t ) ( h ( s , t ) + i ∗ h s + j ∗ h t ) ) ) + c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ( d h ( s , t ) + i ∗ d h s + j ∗ d h t ) = c k ⋅ I k ( i , j ) + c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ( d h ( s , t ) + i ∗ d h s + j ∗ d h t ) \begin{aligned} \boldsymbol{I}_{R}(i, j) \approx \boldsymbol{c}_{k} \cdot & \boldsymbol{I}_{k}\left(\boldsymbol{P}_{k}\left(\boldsymbol{O}_{R}+\overrightarrow{\boldsymbol{r}}_{R}(s, t)\left(h(s, t)+i * h_{s}+j * h_{t}\right)\right)\right)+\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\left(d h(s, t)+i * d h_{s}+j * d h_{t}\right) \\ =& \boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)+\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\left(d h(s, t)+i * d h_{s}+j * d h_{t}\right) \end{aligned} IR(i,j)ck=Ik(Pk(OR+r R(s,t)(h(s,t)+ihs+jht)))+ckh(s,t)Ik(i,j)(dh(s,t)+idhs+jdht)ckIk(i,j)+ckh(s,t)Ik(i,j)(dh(s,t)+idhs+jdht)
    由此,我们可以得到能量函数E的近似表达式为:
    E = ∑ i j k ( I R ( i , j ) − c k ⋅ I k ( i , j ) + c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ( d h ( s , t ) + i ∗ d h s + j ∗ d h t ) ) 2 E=\sum_{i j k}\left(\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)+\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\left(d h(s, t)+i * d h_{s}+j * d h_{t}\right)\right)^{2} E=ijk(IR(i,j)ckIk(i,j)+ckh(s,t)Ik(i,j)(dh(s,t)+idhs+jdht))2

    令:
    b i j k = I R ( i , j ) − c k ⋅ I k ( i , j )    ∇ x = [ d h ( s , t ) , d h s , d h t ] T A i j k = ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) [ 1 i j ] T \boldsymbol{b}_{i j k}=\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)\\\;\\ \nabla \boldsymbol{x}=\left[d h(s, t), d h_{s}, d h_{t}\right]^{T}\\ \boldsymbol{A}_{i j k}=\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]^{T} bijk=IR(i,j)ckIk(i,j)x=[dh(s,t),dhs,dht]TAijk=(ckh(s,t)Ik(i,j))1ijT

    则上式可以表达为:
    E = ∑ i j k ( A i j k ∇ x + b i j k ) 2 E=\sum_{i j k}\left(A_{i j k} \nabla x+b_{i j k}\right)^{2} E=ijk(Aijkx+bijk)2

    E E E关于 ∇ x \nabla x x的梯度,可以得到:
    ∂ E ∂ ∇ x = ∑ i j k A i j k T ( A i j k x + b i j k ) \frac{\partial E}{\partial \nabla \boldsymbol{x}}=\sum_{i j k} \boldsymbol{A}_{i j k}^{T}\left(\boldsymbol{A}_{i j k} \boldsymbol{x}+\boldsymbol{b}_{i j k}\right) xE=ijkAijkT(Aijkx+bijk)

    ∂ E ∂ ∇ x = 0 \frac{\partial E}{\partial \nabla x}=0 xE=0得:
    ( ∑ i j k A i j k T A i j k ) ∇ x + ∑ i j k A i j k T b i j k ) = 0 \left.\left(\sum_{i j k} A_{i j k}^{T} A_{i j k}\right) \nabla x+\sum_{i j k} A_{i j k}^{T} b_{i j k}\right)=0 ijkAijkTAijkx+ijkAijkTbijk=0

    A = ∑ i j k A i j k T A i j k , b = ∑ i j k A i j k b i j k A=\sum_{ijk}A^{T}_{ijk}A_{ijk}, b=\sum_{ijk}A_{ijk}b_{ijk} A=ijkAijkTAijk,b=ijkAijkbijk,则有:
    ∇ x = A − 1 b \nabla x=A^{-1}b x=A1b

    补充:
    A i j k = ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) [ 1 i j ] T    A i j k T A i j k = [ 1 i j ] ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) T ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) [ 1 i j ] T = ∥ c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ∥ 2 [ 1 , i , j i , i 2 , i j j , i j , j 2 ]    A i j k T b i j k = [ 1 i j ] ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) T ( I R ( i , j ) − c k ⋅ I k ( i , j ) ) = s i j k [ 1 i j ] \begin{array}{c} \boldsymbol{A}_{i j k}=\left(\boldsymbol{c}_{k} \cdot \frac{\partial I_{k}(i, j)}{\partial h(s, t)}\right)\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]^{T}\\\;\\ \boldsymbol{A}_{i j k}^{T} \boldsymbol{A}_{i j k}=\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]^{T}=\left\|\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right\|^{2}\left[\begin{array}{ccc} 1, & i, & j \\ i, & i^{2}, & i j \\ j, & i j, & j^{2} \end{array}\right]\\\;\\ \boldsymbol{A}_{i j k}^{T} \boldsymbol{b}_{i j k}=\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right]\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)\right)=s_{i j k}\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right] \end{array} Aijk=(ckh(s,t)Ik(i,j))1ijTAijkTAijk=1ij(ckh(s,t)Ik(i,j))T(ckh(s,t)Ik(i,j))1ijT=ckh(s,t)Ik(i,j)21,i,j,i,i2,ij,jijj2AijkTbijk=1ij(ckh(s,t)Ik(i,j))T(IR(i,j)ckIk(i,j))=sijk1ij

    其中 s i j k = ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) T ( I R ( i , j ) − c k ⋅ I k ( i , j ) ) s_{i j k}=\left(c_{k} \cdot \frac{\partial I_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(I_{R}(i, j)-c_{k} \cdot I_{k}(i, j)\right) sijk=(ckh(s,t)Ik(i,j))T(IR(i,j)ckIk(i,j))是标量,因此可以得到:
    A = ∑ i j k A i j k T A i j k = ∑ i j k ∥ c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ∥ 2 [ 1 , i , j i , i 2 , i j j , i j , j 2 ] b = ∑ i j k A i j k T b i j k = ∑ i j k ( c k ⋅ ∂ I k ( i , j ) ∂ h ( s , t ) ) T ( I R ( i , j ) − c k ⋅ I k ( i , j ) ) [ 1 i j ] \begin{array}{l} \boldsymbol{A}=\sum\limits_{i j k} \boldsymbol{A}_{i j k}^{T} \boldsymbol{A}_{i j k}=\sum\limits_{i j k}\left\|\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right\|^{2}\left[\begin{array}{ccc} 1, & i, & j \\ i, & i^{2}, & i j \\ j, & i j, & j^{2} \end{array}\right]\\ \boldsymbol{b}=\sum \limits_{i j k} \boldsymbol{A}_{i j k}^{T} \boldsymbol{b}_{i j k}=\sum\limits_{i j k}\left(\boldsymbol{c}_{k} \cdot \frac{\partial \boldsymbol{I}_{k}(i, j)}{\partial h(s, t)}\right)^{T}\left(\boldsymbol{I}_{R}(i, j)-\boldsymbol{c}_{k} \cdot \boldsymbol{I}_{k}(i, j)\right)\left[\begin{array}{l} 1 \\ i \\ j \end{array}\right] \end{array} A=ijkAijkTAijk=ijkckh(s,t)Ik(i,j)21,i,j,i,i2,ij,jijj2b=ijkAijkTbijk=ijk(ckh(s,t)Ik(i,j))T(IR(i,j)ckIk(i,j))1ij

    算法流程:

    1. 只进行 h ( s , t ) h(s,t) h(s,t)的优化(迭代4次)
    2. While(iteration <20)
    3.   每间隔5次迭代,进行 h ( s , t ) , h s , h t h(s,t), h_{s}, h_{t} h(s,t)hs,ht的优化,否则仅仅优化深度
    4.   优化完成后跟踪判断每个视角的置信度
    5.   如果视角的置信度太小,或者超过迭代14次尚未收敛,则从局部视角中移除该视角,并重新进行视角选择

    展开全文
    weixin_40957452 2021-09-19 17:42:19
  • 3星
    20.08MB yiyi7392 2018-05-02 17:10:45
  • 3星
    5KB weixin_42233401 2020-10-21 20:20:16
  •  三角测量:已知相机参数和同一个三维点对应多个视角像平面上的匹配点,恢复三维点的坐标(目标是求X的坐标) 已知第i个相机的投影矩阵,其中KiK_{i}Ki​为第 i 个相机的内参矩阵, [pi,ti][p_{i}, t_{i}][pi​,ti​]为...

    SfM(Structure from motion) 是一种传统的三维重建的方法,基于多个视角的图像进行3D重建。从时间系列的多幅2D图像中推算3D信息。

    稀疏点云重建包含以下步骤:

    1. 特征提取与匹配
    2. 获取相机的初始参数
    3. 求解相机姿态
    4. 三角测量得到稀疏点云
    5. 捆绑调整,得到精确的相机参数与稀疏点云坐标

    三角测量

    直接线性变换法

    三角测量:已知相机内外参数和同一个三维点对应多个视角像平面上的同名点坐标,恢复三维点的坐标

    Alt

    已知第i个相机的投影矩阵为 P i P_{i} Pi, 其中 K i K_{i} Ki为第 i 个相机的内参矩阵, [ p i , t i ] [p_{i}, t_{i}] [pi,ti]为外参矩阵,第i个相机的投影方程如(1.1)所示:
    P i = K i [ P i , t i ] = { p i 1 p i 2 p i 3 } (1.1) P_{i}=K_{i}[P_{i}, t_{i}]=\left\{\begin{matrix} p_{i1}\\p_{i2}\\p_{i3} \end{matrix}\right\}\tag{1.1} Pi=Ki[Pi,ti]=pi1pi2pi3(1.1)

    三维点坐标: X i = [ x , y , z , 1 ] T X_{i}=[x,y,z,1]^{T} Xi=[x,y,z,1]T, 在第i个视角的投影的图像(归一化像平面)坐标为: x i ∗ = [ x i , y i , 1 ] T x_{i}^{*}=[x_{i},y_{i},1]^{T} xi=[xi,yi,1]T
    根据投影方程: x i ∗ = P i X (1.2) x^{*}_{i}=P_{i}X\tag{1.2} xi=PiX(1.2)

    (1.2)两侧同时左叉乘 x i ∗ x^{*}_{i} xi得: 0 = x i ∗ × ( P i X ) = x i ∗ × ( { p i 1 p i 2 p i 3 } X ) = { x i y i 1 } × { p i 1 X p i 2 X p i 3 X } (1.3) 0=x^{*}_{i}\times(P_{i}X)=x^{*}_{i}\times(\left\{\begin{matrix} p_{i1}\\p_{i2}\\p_{i3} \end{matrix}\right\}X)=\left\{\begin{matrix}x_{i}\\y_{i}\\1\end{matrix}\right\}\times \left\{\begin{matrix} p_{i1}X\\p_{i2}X\\p_{i3}X \end{matrix}\right\}\tag{1.3} 0=xi×(PiX)=xi×(pi1pi2pi3X)=xiyi1×pi1Xpi2Xpi3X(1.3)

    即:
    x i ( p i 3 X ) − P i 1 X = 0 y i ( p i 3 X ) − P i 2 X = 0 x i ( p i 2 X ) − y 1 ( P i 1 X ) = 0 (1.4) x_{i}(p_{i3}X)-P_{i1}X=0\\y_{i}(p_{i3}X)-P_{i2}X=0\\x_{i}(p_{i2}X)-y_{1}(P_{i1}X)=0\tag{1.4} xi(pi3X)Pi1X=0yi(pi3X)Pi2X=0xi(pi2X)y1(Pi1X)=0(1.4)

    由(1.4)得, 第三个等式与前两个等式线性相关,将其写成矩阵形式如(1.5)所示:
    { x i P i 3 − P i 1 y i P i 3 − P i 2 } X = 0 (1.5) \left\{\begin{matrix} x_{i}P_{i3}-P{i1}\\y_{i}P_{i3}-P{i2} \end{matrix}\right\}X=0\tag{1.5} {xiPi3Pi1yiPi3Pi2}X=0(1.5)

      已知 X = [ x i , y i , z i , 1 ] T X=[x_{i},y_{i},z_{i},1]^{T} X=[xi,yi,zi,1]T, 显然X的自由度为3,由(1.5)可知X点投影的一个视角可以提供两个约束,至少需要两个以上视角才能求解X的三维坐标. 如果同一个三维点X,有N个视角的图像的投影坐标 [ ( x 1 , y 1 ) , . . . ( x N , y N ) ] [(x_{1},y_{1}),...(x_{N},y_{N})] [(x1,y1),...(xN,yN)]与对应N个相机内外参数,通常使用最小二乘的方式来求解,或者采用Ransac结合最小二乘法求解X的坐标. A X = 0 A = [ x 1 P 13 − P 11 y 1 P 13 − P 12 . . . x i P i 3 − P i 1 y i P i 3 − P i 2 . . . x N P N 3 − P N 1 y N P N 3 − P N 2 ] , N ≥ 2 (1.6) AX=0\\A=\begin{bmatrix}x_{1}P_{13}-P_{11}\\y_{1}P_{13}-P_{12}\\...\\x_{i}P_{i3}-P_{i1}\\y_{i}P_{i3}-P_{i2}\\...\\x_{N}P_{N3}-P_{N1}\\y_{N}P_{N3}-P_{N2}\end{bmatrix}, N\ge2 \tag{1.6} AX=0A=x1P13P11y1P13P