解决：如何在opencv中得到与matlab立体标定一样的矫正图？（python版opencv）

目的：采用一样的标定参数，matlab中和opencv中的立体矫正图像是一样的吗？不一样的话怎么让它们一样？
结论：不一样。后文为解决方案。
原因：注意matlab的标定结果在matlab中的用法和在opencv中的用法不一样，主要原因是matlab的rectifyStereoImages函数和opencv的cv2.stereoRectify函数的计算结果不一样导致的。照这个思路，把matlab的rectifyStereoImages的结果导入opencv、而不用opencv的cv2.stereoRectify就可以了。
不想看对比细节可直接看3. 好像发现了一些线索
软件版本：Matlab2024a + OpenCV4.9.0

问题的提出
在matlab的stereo camera calibrator中对棋盘格标定后，可以点击“Show Rectified”一键得到矫正图
在这里插入图片描述
但是把matlab的参数导入到opencv中后，怎么得到和matlab一样的结果呢？

1. 用标定后的参数在matlab中直接显示立体矫正后的左右目图像

在matlab的stereo camera calibrator中对棋盘格标定后，有两种结果，一种直接在stereoParams中查看，另一种用函数转换成opencv的格式，这里两种都展示：

1.1 在stereoParams中直接查看

% stereoParams.CameraParameters1.intrinsics中查看左相机
K =[5701.4907,0,1387.8609;0,5705.7330,1050.1707;0,0,1]
RadialDistortion = [-0.0621,0.2954,2.6695]
TangentialDistortion =[0.0010,-0.0004]
%写成[k1, k2, p1, p2, k3]格式
distCoeffs = [-0.0621,0.2954,0.0010,-0.0004,2.6695]
% stereoParams.CameraParameters2.intrinsics中查看右相机
K =[5596.4954,0,1062.2055;0,5620.0099,1016.7026;0,0,1]
RadialDistortion = [-0.0358,-0.0958,9.2830]
TangentialDistortion =[0.0001,0.0034]
%写成[k1, k2, p1, p2, k3]格式
distCoeffs = [-0.0358,-0.0958,0.0001,0.0034,9.2830]
% 以及旋转矩阵、平移矩阵（这两个是完全一样的）
stereoParams.PoseCamera2.R =[0.9500,-0.0045,0.3121;0.0035,1.0000,0.0038;-0.3121,-0.0025,0.9500]
stereoParams.PoseCamera2.Translation =[  -35.0696   -0.0101    5.1674]

1.2 在matlab中转换成opencv格式

函数：stereoParametersToOpenCV

[intrinsicMatrix1,distortionCoefficients1,...intrinsicMatrix2,distortionCoefficients2,...rotationOfCamera2,translationOfCamera2] ...= stereoParametersToOpenCV(stereoParams)
% 输出
intrinsicMatrix1 =[5701.4907,0,1386.8609;0,5705.7330,1049.1707;0,0,1]
distortionCoefficients1 = [-0.0621    0.2954    0.0010   -0.0004    2.6695]
intrinsicMatrix2 =[5596.4954,0,1061.2055;0,5620.0099,1015.7026;0,0,1]
distortionCoefficients2 =[-0.0358   -0.0958   0.0001    0.0034   9.2830]
rotationOfCamera2 =[0.9500   -0.0045	0.31210.0035    1.0000	0.0038-0.3121   -0.0025	0.9500]
translationOfCamera2 =[  -35.0696   -0.0101    5.1674]

1.3 两种方法的差异

只有内参矩阵的主点坐标不一样，畸变参数、旋转矩阵、平移向量都是一样的~
差值为1，原因是主点-1了，如下：cx = intrinsics.PrincipalPoint(1) - 1;

%在函数function [intrinsicMatrix, distortionCoefficients] = cameraIntrinsicsToOpenCV(intrinsics)中
function distortionCoefficients = getOCVDistortionCoefficients(intrinsics)    % Distortion Coefficients in OpenCV [k1 k2 p1 p2 k3]distortionCoefficients = zeros(1,5);if length(intrinsics.RadialDistortion) == 3distortionCoefficients([1,2,5]) = intrinsics.RadialDistortion;elsedistortionCoefficients([1,2]) = intrinsics.RadialDistortion;enddistortionCoefficients([3,4])   = intrinsics.TangentialDistortion;
endfunction intrinsicMatrix = getOCVIntrinsicMatrix(intrinsics)% Focal length.fx = intrinsics.FocalLength(1);fy = intrinsics.FocalLength(2);% Principal point.cx = intrinsics.PrincipalPoint(1) - 1;cy = intrinsics.PrincipalPoint(2) - 1;        % Construct OpenCV's intrinsic matrix.intrinsicMatrix = [fx   0  cx;0  fy  cy;0   0  1];
end

1.4 在matlab中用stereoParams映射矫正后的图像

close all;  clc;I1 = imread('D:\StereoRectify\Calibration\left\0001.jpg');%读取左右图片
I2 = imread('D:\StereoRectify\Calibration\right\0001.jpg');[J1, J2] = rectifyStereoImages(I1, I2, stereoParams);figure; subplot(2,1,1); imshow(J1);
subplot(2,1,2); imshow(J2);

这里J1，J2的图像尺寸都是4026×2128。
如果在opencv里面矫正后的图像尺寸也是这个值，说明两个方法取得一样的效果；反之，就有问题！
当然，这个验证方法过于简单粗暴，实际上把两种方法的图像做减法也能对比的~~
在这里插入图片描述

但是把matlab的参数导入到opencv中后，怎么得到和matlab一样的结果呢？

2. 把matlab标定的参数导入python的opencv，然后计算立体矫正后的左右目图像

这段可以不看~~怎么解释validPixROI这个参数呢？
在这里插入图片描述
网上有这么说的：

validPixROI1：一个最多地包含有效像素的长方形。（左目图像）
validPixROI2：一个最多地包含有效像素的长方形。（右目图像）

这个参数的第一个值是有效图像最左侧非零点的横坐标；
后面三个参数我无法解释，对着图像来看上/下/左/右/横/纵坐标都不是-_-|||。
opencv官方文档也没写，我自己写一个计算图像的非零区域的代码吧！

跳到这里：这里是正文~~

2.1 代码：

import camera_config
import cv2
import matplotlib.pyplot as plt
from PIL import Image
import numpy as npsize = camera_config.size  # (8000, 3000) # 图像尺寸
left_camera_matrix = camera_config.left_camera_matrix
left_distortion = camera_config.left_distortion
right_camera_matrix = camera_config.right_camera_matrix
right_distortion = camera_config.right_distortion
R = camera_config.R
T = camera_config.T# 进行立体更正
R1, R2, P1, P2, Q, validPixROI1, validPixROI2 = cv2.stereoRectify(left_camera_matrix, left_distortion, right_camera_matrix, right_distortion, size, R, T)
# 输出参数：
print('***** R1 *****\n',R1)    # R1：矫正旋转矩阵。将第一个相机坐标系下未矫正的点变换到第一个相机矫正坐标系下，即 R_{左矫正坐标系}{左未矫正坐标系}
print('***** P1 *****\n',P1)    # P1：3x4左相机投影矩阵。将左矫正坐标系下的点投影到左矫正坐标系图像平面坐标系。
print('***** R2 *****\n',R2)    # R2：矫正旋转矩阵。将第二个相机坐标系下未矫正的点变换到第二个相机矫正坐标系下，即 R_{右矫正坐标系}{右未矫正坐标系}
print('***** P2 *****\n',P2)    # P2：3x4右相机投影矩阵。将左矫正坐标系下的点投影到右矫正坐标系图像平面坐标系。
print('***** Q *****\n',Q)      # Q：4x4的视差深度映射矩阵。
# 计算更正map
left_map1, left_map2 = cv2.initUndistortRectifyMap(left_camera_matrix, left_distortion, R1, P1, size, cv2.CV_16SC2)
right_map1, right_map2 = cv2.initUndistortRectifyMap(right_camera_matrix, right_distortion, R2, P2, size, cv2.CV_16SC2)
img_l = np.array(Image.open(r'Calibration/left/0001.jpg')) 
img_r = np.array(Image.open(r'Calibration/right/0001.jpg')) 
# 根据更正map对图片进行重构
img1_rectified = cv2.remap(img_l, left_map1, left_map2, cv2.INTER_LINEAR)
img2_rectified = cv2.remap(img_r, right_map1, right_map2, cv2.INTER_LINEAR)
# 转成灰度图，求非零图像区域，然后求图像区域的上/下/左/右/横/纵坐标
img1_rectified_gray = cv2.cvtColor(img1_rectified, cv2.COLOR_RGB2GRAY)
img2_rectified_gray = cv2.cvtColor(img2_rectified, cv2.COLOR_RGB2GRAY)
idxNonZero=cv2.findNonZero(img1_rectified_gray).squeeze(1)
l1, r1, t1, b1 = np.min(idxNonZero[:,0]), np.max(idxNonZero[:,0]), np.min(idxNonZero[:,1]), np.max(idxNonZero[:,1])
idxNonZero=cv2.findNonZero(img2_rectified_gray).squeeze(1)
l2, r2, t2, b2 = np.min(idxNonZero[:,0]), np.max(idxNonZero[:,0]), np.min(idxNonZero[:,1]), np.max(idxNonZero[:,1])
plt.figure(1)
plt.subplot(211)
plt.imshow(img1_rectified)
plt.subplot(212)
plt.imshow(img2_rectified)
plt.figure(2)
plt.subplot(211)
plt.imshow(img1_rectified[min(t1,t2):max(b1,b2), l2:r1])
plt.subplot(212)
plt.imshow(img2_rectified[min(t1,t2):max(b1,b2), l2:r1])
# 裁剪后的图像尺寸
print(f'image size: {r1-l2}x{max(b1,b2)-min(t1,t2)}')plt.show()

结果：

***** R1 *****[[ 0.9853897  -0.00381456  0.1702721 ][ 0.00352377  0.99999177  0.00200993][-0.17027837 -0.00138056  0.98539503]]
***** P1 *****[[5662.8715        0.         4312.4698143     0.        ][   0.         5662.8715     1261.70807171    0.        ][   0.            0.            1.            0.        ]]
***** R2 *****[[ 0.98931806  0.00028444 -0.14577278][-0.00053314  0.99999847 -0.00166697][ 0.14577209  0.00172688  0.98931669]]
***** P2 *****[[   5662.8715           0.            4312.4698143  -200738.92039502][      0.            5662.8715        1261.70807171       0.        ][      0.               0.               1.               0.        ]]
***** Q *****[[    1.             0.             0.         -4312.4698143 ][    0.             1.             0.         -1261.70807171][    0.             0.             0.          5662.8715    ][    0.             0.             0.02821013    -0.        ]]
image size: 4074x2152

2.2 参数解释

参见opencv官方文档中函数stereoRectify()的解释

P1, P2的含义及其用法

P1：3x4左相机投影矩阵。将左矫正坐标系下的三维点投影到左矫正坐标系图像平面坐标系的像素点
$}\left[ {\begin{array}{cc} f&0&{{c_{x1}}}&0\\ 0&f&{{c_y}}&0\\ 0&0&1&0 \end{array}} \right]$

P2：3x4右相机投影矩阵。将左矫正坐标系下的三维点投影到右矫正坐标系图像平面坐标系的像素点。 $T_x$ 是两个相机之间的水平偏移值。我这里算了一下，实际上就是matlab标定的平移向量[ -35.0696 -0.0101 5.1674] 的长度=35.4483的负数，用P2的两个参数 ${T_x} \cdot f$ =-200738.92039502； $f$ =5662.8715； -200738.92039502/5662.8715 = -35.4483 可以验证。
$}\left[ {\begin{array}{cc} f&0&{{c_{x2}}}&{{T_x} \cdot f}\\ 0&f&{{c_y}}&0\\ 0&0&1&0 \end{array}} \right]$

对于三维物理空间点 $(X, Y, Z)$ ，都可用P1和P2，计算该点对应的图像像素坐标 $(x, y)$ ，（ $w$ 是尺度因子）
$w\left[ {\begin{array}{cc} x\\ y\\ 1 \end{array}} \right] = P\left[ {\begin{array}{cc} X\\ Y\\ Z\\ 1 \end{array}} \right]$

Q的含义及其用法

Q：4x4的视差深度映射矩阵。
$}\left[ {\begin{array}{cc} 1&0&0&{ - {c_{x1}}}\\ 0&1&0&{ - c{}_y}\\ 0&0&0&f\\ 0&0&{ - \frac{1}{{{T_x}}}}&{\frac{{{c_{x1}} - {c_{x2}}}}{{{T_x}}}} \end{array}} \right]$
为什么把Q叫视差深度映射矩阵呢？因为采用下式，Q能够将单通道视差图转换为表示 3D 表面的 3 通道图（就是x,y,z坐标值）。对于每个像素坐标 $(x, y)$ 及其视差 $d = d i s p a r i t y (x, y)$ ，都能计算三维物理空间中的一个点 $(X, Y, Z)$ ：
$\left[ {\begin{array}{c} X\\ Y\\ Z\\ W \end{array}} \right] = Q\left[ {\begin{array}{c} x\\ y\\ disparity(x,y)\\ 1 \end{array}} \right]$
上式中的W好像没有什么实际意义。。。

2.3 图像结果

Figure1是8000×3000尺寸的图像
在这里插入图片描述
Figure2是裁掉周围黑框的图像

opencv的image size: 4074x2152为什么和matlab的图像尺寸4026×2128不一样啊？！

我裂开！！！

哪位小伙伴知道为什么？欢迎讨论！！！
直观的对比opencv和matlab的结果
在这里插入图片描述

3. 好像发现了一些线索 ~_~

matlab自带的rectifyStereoImages(I1, I2, stereoParams)的可选输出重投影矩阵reprojectionMatrix，和opencv的Q矩阵是同样的含义，但是，结果却不一样！！！

[J1, J2, reprojectionMatrix, camMatrix1, camMatrix2, R1, R2] = rectifyStereoImages(I1, I2, stereoParams);

reprojectionMatrix = [1	0	0	-1943.86090	1	0	-1066.43670	0	0	5596.49540	0	0.0282	0]
%   Notes
%   -----
%   - reprojectionMatrix is represented as a 4-by-4 matrix:
%       [1    0   0     -cx
%        0    1   0     -cy
%        0    0   0     f
%        0    0   1/b   0],
%     where f and [cx, cy] are the focal length and principal point of
%     rectified camera 1, respectively. b is the baseline of the virtual 
%     rectified stereo camera.
camMatrix1 =[5596.4954	0	1943.8609	-198385.9933;0	5596.4954	1066.4367	0;0	0	1	0]
camMatrix2 =[1	0	0	0;0	1	0	0;0	0	1	0]
% R1、R2和前面的一样
R1 =[0.9854	-0.0038	0.17030.0035	1.0000	0.0020-0.1703	-0.0014	0.9854]
R2 =[0.9893	0.0003	-0.1458-0.0005	1.0000	-0.00170.1458	0.0017	0.9893]
%   - Use camMatrix1 and camMatrix2 to project 3-D world points in the 
%     rectified camera 1's coordinate system into the image plane of J1 and 
%     J2, respectively. R1 and R2 bring 3-D points in the unrectified camera's 
%     coordinate system to points in the rectified camera's coordinate system 
%     for camera 1 and camera 2, respectively.

在python中对前面的2.1 代码作出如下修改就可以了：

# 原来的
# 进行立体更正
R1, R2, P1, P2, Q, validPixROI1, validPixROI2 = cv2.stereoRectify(left_camera_matrix, left_distortion, right_camera_matrix, right_distortion,size, R, T)# 把上面这句删除，改成matlab的rectifyStereoImages输出值
Q = reprojectionMatrix
P1 = camMatrix1
P2 = camMatrix2
R1 = R1
R2 = R2

这样opencv输出的image size: 4025x2111，和matlab的图像尺寸4026×2128相当接近了！！！直观对比如下 ↓↓
在这里插入图片描述

但是，左图是可以的，右图是黑色的，怎么办啊？？！！

尝试：按照前面P2的定义，对P2重新赋值

P2 = P1
P2[0,3] = np.linalg.norm(T)*P1[0,0]
# 输出P2
***** P2 *****[[  5596.495           0.           1943.861      198385.9892][     0.           5596.495        1066.437           0.        ][     0.              0.              1.              0.        ]]image size: 4025x2127