<think>我们有两张影像（tif格式）和对应的投影矩阵。投影矩阵文件提供了两个矩阵，每个矩阵有4行，每行有4个数值（实际上我们看到的是3行，但每行有4个数值，所以是3x4的矩阵，通常用于表示齐次坐标下的3x4投影矩阵，即P = [M | t]，其中M是3x3，t是3x1）。
根据核线影像生成的原理，我们知道核线影像是为了消除立体像对中的y方向视差，使得同名点位于同一行（或同一列）上。通常，核线影像分为水平核线影像（核线平行于x轴，即同名点在同一行）和竖直核线影像（核线平行于y轴，即同名点在同一列）。而最接近原始图像的核线影像通常是指在重采样过程中失真最小的核线影像。
流程：
1. 读取两张影像和对应的投影矩阵（P1, P2）。
2. 根据投影矩阵计算基础矩阵F（Fundamental Matrix）。基础矩阵可以通过两个投影矩阵计算：F = [e2]x * P2 * P1^+，其中P1^+是P1的伪逆，e2是第二个相机的光心在第一个相机图像上的投影（即极点），[e2]x是e2的叉乘矩阵。
但是更常见的方法是：已知两个投影矩阵P1和P2，则F = [P2]_{3x3} * [P1]_{3x3}^{-T} * [e2]x，但实际上有更直接的方式：我们可以通过以下步骤计算：
- 计算两个相机的光心：C1 = -M1^{-1} * t1, C2 = -M2^{-1} * t2，其中P1 = [M1 | t1]，P2 = [M2 | t2]（注意这里的M1是3x3，t1是3x1）
- 计算第二个相机的光心C2在第一个相机图像上的投影：e1 = P1 * [C2; 1]（注意齐次坐标，C2是3维，加1变成4维）
- 计算基础矩阵F = [e1]x * P1 * P2^+，其中P2^+是P2的伪逆（因为P2是3x4，伪逆是4x3的矩阵）。
然而，另一种标准方法是：F = [e1]x * (P1 * P2^+)，其中P2^+是P2的伪逆（Moore-Penrose伪逆）。
但是，由于我们已知投影矩阵，更简单的方法是使用对极几何关系。实际上，基础矩阵F满足：x2^T * F * x1 = 0，其中x1和x2是两幅图像上的对应点。
3. 根据基础矩阵F，我们可以计算核线变换矩阵（H1和H2），使得变换后的图像对满足核线条件。常用的方法有：
- 使用Hartley的归一化8点法计算基础矩阵（但这里我们已经有了投影矩阵，所以可以直接计算F，不需要匹配点）。
- 或者，我们可以通过计算两个图像的核线几何，然后利用单应矩阵将核线映射到水平（或竖直）方向。
4. 具体地，对于水平核线影像，我们希望将核线映射为水平线。这可以通过对每幅图像应用一个单应变换来实现。常用的方法是：
- 选择第一幅图像的一个变换H1（例如，我们可以选择H1使得第一幅图像中的核线变为水平，同时保持图像尽可能不变形）。
- 然后，第二幅图像的变换H2可以通过H2 = H1 * inv(F) 或者通过其他约束得到？实际上，Hartley提出了一种方法：首先计算第一幅图像的一个单应H1（例如，将极点映射到无穷远点[1,0,0]^T，这样核线就变成水平的了），然后第二幅图像的单应H2可以通过H2 = [e2]x * F + e2 * v^T（其中v是任意向量）来构造，但这样比较复杂。
5. 实际上，OpenCV提供了函数`cv::stereoRectify`来直接计算两个图像的校正变换（即单应矩阵），该函数需要输入两个相机的内参、外参（旋转和平移）等。但是，我们这里只有投影矩阵P1和P2，所以需要从投影矩阵中分解出内参和旋转平移。
6. 另一种思路：由于我们已知投影矩阵，我们可以分解投影矩阵得到内参矩阵K和旋转矩阵R以及平移向量t（通过RQ分解）。然后，使用这些参数进行立体校正。
7. 步骤：
a. 读取两个投影矩阵P1和P2（3x4）。
b. 分解投影矩阵：将投影矩阵分解为K1, R1, t1和K2, R2, t2。分解方法：使用QR分解（或RQ分解）分解投影矩阵的前3x3部分（记为M），因为P = K [R | t]，其中K是上三角矩阵（内参），R是旋转矩阵（正交矩阵），t是平移向量。注意：分解时，需要将M分解为K和R，使得M=K*R，其中K是上三角矩阵，R是正交矩阵。
c. 然后，我们可以使用分解得到的参数（K1, R1, t1和K2, R2, t2）以及两个相机之间的相对位姿（R和t）来调用立体校正函数。实际上，两个相机之间的旋转矩阵R和平移向量t可以通过以下方式得到：
R = R2 * R1^T
t = t2 - R * t1 （注意这里t1和t2是相机在世界坐标系中的位置？实际上，我们需要的是从第一个相机到第二个相机的变换：R12和t12）
更准确地说，如果我们有第一个相机的投影矩阵P1 = K1 [I | 0]（以第一个相机为世界坐标系），而第二个相机的投影矩阵P2 = K2 [R | t]（其中R和t是第一个相机坐标系到第二个相机坐标系的变换），那么我们可以这样分解：
首先，将P1分解：P1 = K1 [I | 0]
然后，P2 = K2 [R | t] = K2 [R | -R*C2]，其中C2是第二个相机在第一个相机坐标系中的位置（即光心）。
但是，我们现有的投影矩阵P1和P2是相对于某个世界坐标系的，所以我们需要先定义世界坐标系。通常，我们可以将第一个相机的位置作为世界坐标系原点，然后分解P1得到K1和[R1|t1]（此时R1=I，t1=0），但这里P1和P2都是相对于同一个世界坐标系的。
因此，我们首先分解P1：P1 = [M1 | t1] = [K1 * R1 | t1']，其中t1' = -K1*R1*C1，而C1是第一个相机的光心（在世界坐标系下）。同样分解P2：P2 = [K2 * R2 | t2']。
然后，两个相机之间的相对旋转为：R = R2 * R1^T，相对平移为：t = t2' - R * t1'（注意这里的平移向量是在第一个相机坐标系下的？）实际上，我们需要的是从第一个相机坐标系到第二个相机坐标系的变换，即点P在第一个相机坐标系下的坐标X1，在第二个相机坐标系下为X2 = R * X1 + t。
然而，在立体校正中，我们通常需要的是两个相机之间的相对旋转（R1_to2）和平移（t1_to2）。我们可以这样计算：
R1_to2 = R2 * R1^T
t1_to2 = t2' - R1_to2 * t1' ? 这个关系并不直接。
实际上，更常用的方法是：我们已知两个投影矩阵P1和P2，那么基础矩阵F可以通过以下公式计算：
F = [e2]_x * P2 * P1^+
其中，e2是第二个相机的光心在第一个相机图像上的投影（即极点），P1^+是P1的伪逆。
但为了使用OpenCV的立体校正函数，我们更倾向于分解出内参和相对位姿。
8. 由于分解投影矩阵比较复杂，而且容易出错，我们可以考虑使用另一种方法：直接使用基础矩阵F和两个极点来构造核线变换矩阵。Hartley在论文《Theory and Practice of Projective Rectification》中提出了一种方法，称为“无标定立体校正”，只需要基础矩阵F和一个匹配点（或极点）即可。
9. 具体步骤如下（Hartley方法）：
a. 计算基础矩阵F（使用投影矩阵）。
b. 计算两个极点：e1和e2（e1是第二个相机在第一个相机图像上的极点，e2是第一个相机在第二个相机图像上的极点）。可以通过F*e1=0和F^T*e2=0计算，其中e1是F的右零向量（即F^T的零空间），e2是F的左零向量（即F的零空间）。
c. 对于第一幅图像，定义变换H1，将极点e1映射到无穷远点[1,0,0]^T。这样，变换后的核线将变成水平线。H1可以构造为：
H1 = [I | 0] * G （这里不是直接这样，而是构造一个3x3矩阵）
具体构造方法：
- 首先，将图像坐标系的原点平移到e1（这样极点就在原点）：T = [1,0,-e1_x; 0,1,-e1_y; 0,0,1]
- 然后，计算一个旋转矩阵R，使得平移后的极点方向（即T*e1）旋转到x轴上。但注意，平移后的极点是(0,0,1)（齐次坐标），所以我们需要旋转？实际上，平移后极点在原点，但我们需要将核线旋转成水平。然而，Hartley提出了一种更简单的方法：构造一个矩阵将极点映射到[1,0,0]^T，同时最小化图像变形。
Hartley的构造方法：
H1 = [a, b, c;
d, e, f;
g, h, i]
要求H1 * e1 = [1,0,0]^T（齐次坐标下，即H1*e1与[1,0,0]^T共线）。因此，我们可以任意选择H1，只要满足H1*e1=[1,0,0]^T（乘以一个常数因子）。同时，为了最小化图像变形，我们要求H1尽可能接近仿射变换（即最后一行[g,h,i]尽可能接近[0,0,1]），并且对图像进行平移和旋转，使得图像内容在变换后保持原状。
实际上，Hartley提出了一种具体的构造方法（见论文）：
H1 = [1, 0, 0;
0, 1, 0;
-1/e1_x, -1/e1_y, 1] # 这个不对
更准确地说，我们构造一个将e1映射到[f,0,1]^T（然后归一化为[f,0,1]^T，再乘以一个因子使其变成[1,0,0]^T）的矩阵。但这里我们采用以下步骤：
Step 1: 构造一个平移矩阵T，将e1平移到原点：T = [1,0,-e1_x; 0,1,-e1_y; 0,0,1]
Step 2: 构造一个旋转矩阵R，将平移后的点（原点）旋转到x轴上的某个点？实际上，平移后极点位于原点，但核线是过原点的直线。我们希望将核线旋转成水平。然而，Hartley建议使用以下变换：
G = [1,0,0; 0,1,0; -1/f_x, -1/f_y, 1] # 这只是一个例子
实际上，Hartley的方法中，H1由两个矩阵相乘得到：H1 = H_A * H_R
H_A = [1,0,0; 0,1,0; -1/e1_x, -1/e1_y, 1] # 将e1映射到[0,0,1]^T？不对。
另一种常见的构造方法（来自Hartley的论文）：
H2 = [I | 0] # 不对
这里我们采用标准方法（来自OpenCV的文档和Hartley的论文）：
定义H1 = G * T，其中T是平移矩阵（将e1移到原点），G是一个3x3矩阵，定义为：
G = [1, 0, 0;
0, 1, 0;
a, b, c]
要求G将原点映射到[1,0,0]^T（即齐次坐标下[0,0,1]映射到[1,0,0]^T），所以有：c * [a, b, c]^T = [1,0,0]^T -> 显然不可能。因此，我们需要重新考虑。
实际上，Hartley的论文中给出了一个具体的构造步骤（无迭代）：
1. 将e1平移到原点：T = [1,0,-e1_x; 0,1,-e1_y; 0,0,1]
2. 旋转使得e1方向与x轴重合？但e1现在在原点，所以旋转没有意义。因此，我们跳过旋转（或者用旋转来调整图像方向，但为了最小化变形，我们通常不旋转，或者旋转角度很小）。
3. 构造一个矩阵M，使得极点映射到无穷远点[1,0,0]^T：M = [1,0,0; 0,1,0; a,b,1]（注意：这里我们希望M将原点（齐次坐标[0,0,1]）映射到[0,0,1]（因为平移后极点在原点），但我们希望映射后极点在无穷远，即齐次坐标的第三个分量为0。因此，我们需要一个映射使得[0,0,1]映射到[0,0,0]？这不可能。所以我们需要重新考虑：我们希望将e1（在平移后是[0,0,1]）映射到[1,0,0]^T（无穷远点）。因此，我们需要一个矩阵M，使得M*[0,0,1]^T = [0,0,0]^T？这不对，因为齐次坐标乘以矩阵后，[0,0,1]^T乘以M后应该是[ M_{31}, M_{32}, M_{33} ]，我们要求这个向量与[1,0,0]^T共线，即M_{31}:M_{32}:M_{33} = 1:0:0。因此，我们可以取M_{31}=1, M_{32}=0, M_{33}=0。但这样M的最后一行是[1,0,0]，那么M为：
M = [1,0,0;
0,1,0;
1,0,0]
但是，这个矩阵不可逆（秩为2），所以我们需要重新考虑。
正确的构造（Hartley论文中的公式）：
H = [a, b, c;
d, e, f;
g, h, i]
要求H * e1 = [1,0,0]^T（即H将e1映射到无穷远点[1,0,0]^T）。因此，我们取任意一个满足H*e1 = k*[1,0,0]^T（k非零）的H。同时，为了最小化变形，Hartley建议将H分解为仿射部分和扭曲部分，并让仿射部分尽可能保持图像不变。
具体构造（一种简单方法）：
H1 = [e1[0], e1[1], e1[2];
-e1[1], e1[0], 0;
-e1[0]*e1[2], -e1[1]*e1[2], e1[0]^2+e1[1]^2]
然后归一化，使得H1[0]的模为1。但这种方法可能不是最优的。
另一种常见做法（来自网络）：
H1 = [1, 0, 0;
0, 1, 0;
-1/e1_x, -1/e1_y, 1] # 这显然不对，因为当e1_x或e1_y为0时会出现除0错误。
由于时间关系，我们采用OpenCV的`stereoRectify`函数，但需要分解投影矩阵。因此，我们选择分解投影矩阵。
10. 分解投影矩阵（RQ分解）：
将投影矩阵P分解为P = K [R | t]，其中K是上三角矩阵（内参），R是旋转矩阵（正交矩阵），t是平移向量。
分解步骤：
a. 取P的前3列（3x3矩阵）记为M。
b. 对M进行RQ分解（即M = K * R，其中K是上三角矩阵，R是正交矩阵）。注意，RQ分解是QR分解的转置：M^T = Q * R^T，然后K = R^T, R = Q^T？或者直接使用RQ分解函数。
在OpenCV中，我们可以使用`cv::RQDecomp3x3`函数进行RQ分解。
分解后，我们得到K1, R1（对于第一个投影矩阵），同样分解第二个投影矩阵得到K2, R2。
然后，计算两个相机之间的相对旋转和平移：
- 第一个相机的位姿：旋转矩阵R1，平移向量t1（由P1的第四列：t1 = K1^{-1} * p14，其中p14是P1的第四列的前三个元素？因为P1 = K1 [R1 | t1] = K1*R1 | K1*t1，所以t1 = K1^{-1} * p14）
- 第二个相机的位姿：R2, t2（同理：t2 = K2^{-1} * p24）
- 相对旋转：R = R2 * R1^T
- 相对平移：t = t2 - R * t1 （注意：这里t是在第一个相机坐标系下的平移向量）
然后，我们可以使用`cv::stereoRectify`函数计算校正变换矩阵（单应矩阵）H1和H2。
11. 代码步骤：
a. 读取两个tif图像（使用GDAL或OpenCV，但tif可能包含地理信息，所以建议用GDAL，但为了简单，我们假设是普通图像，用OpenCV的imread读取，注意可能有多波段）。
b. 读取投影矩阵文件（两个3x4矩阵）。
c. 分解两个投影矩阵（使用RQ分解）。
d. 计算两个相机之间的相对旋转R和平移向量t。
e. 调用`cv::stereoRectify`函数，得到校正变换矩阵（单应矩阵）H1和H2，以及校正后的投影矩阵P1_new和P2_new（可选）。
f. 对两个图像应用单应变换，得到校正后的图像（水平核线影像）。
g. 竖直核线影像：实际上，立体校正后通常得到水平核线影像。竖直核线影像可以通过旋转90度得到，但这不是通常的做法。通常我们只做水平校正。所以竖直核线影像可能是指将核线变成竖直方向，这可以通过将图像旋转90度后再做校正，然后再旋转回来？或者改变校正方向。但通常我们只做水平校正。
h. 最接近原始图像的核线影像：通常，我们选择第一幅图像的校正变换H1尽可能接近单位矩阵（即最小化变形），或者使用参数设置（在stereoRectify中设置alpha参数，0表示不保留黑边，1表示保留所有像素，通常取-1表示自动选择，但我们可以调整alpha使得图像变形最小）。
12. 关于竖直核线影像：我们可以通过交换图像的x和y坐标（即转置）来模拟竖直核线，但这样会改变图像的方向。另一种方法是：在立体校正时，我们要求校正后的核线是竖直的（即同名点在同一列）。在OpenCV的`stereoRectify`函数中，我们可以通过设置flags参数来指定校正后的核线方向（水平或竖直），但默认是水平的。我们可以通过设置cv::RECTIFY_VERTICAL来获得竖直核线影像（但注意，OpenCV的文档中并没有这个选项，实际上RECTIFY_VERTICAL是用于鱼眼相机的，而且通常我们使用水平核线）。因此，我们可能需要两次校正：一次得到水平核线，一次得到竖直核线（通过交换两个相机的图像角色？或者旋转图像？）。或者，我们可以旋转图像90度，然后进行水平校正，再旋转回来，这样核线就变成竖直的了。
13. 由于时间限制，我们只实现水平核线影像。对于竖直核线影像，我们可以通过将图像旋转90度，然后进行水平校正，最后再旋转-90度来实现。这样，在旋转后的图像上，水平核线就对应原始图像的竖直核线。
14. 最接近原始图像的核线影像：我们可以对第一幅图像使用单位矩阵（即不校正），但这不满足核线条件。所以，我们可以在校正时调整`stereoRectify`中的alpha参数，使得图像裁剪最少（即保留更多有效区域，同时变形较小）。alpha=0表示校正后的图像会进行缩放和平移，使得所有像素都是有效的（即没有黑边），但可能会丢失一些像素；alpha=1表示保留所有原始图像像素（包括黑边）。我们可以设置alpha为0.5，或者尝试不同的值，然后选择看起来变形最小的那个。或者，我们可以计算两个校正变换（H1和H2）的变形程度（例如，计算变换后图像的纵横比变化、角度变化等），然后选择变形最小的那个。但通常，我们使用alpha=-1（自动选择）来平衡有效区域和变形。
15. 代码实现步骤：
a. 读取图像：img1, img2（使用OpenCV的imread，注意路径）。
b. 读取投影矩阵文件，解析出两个3x4矩阵P1和P2（注意：文件中的矩阵是行优先，每行4个数值）。
c. 对P1和P2进行分解（得到K1, R1, t1和K2, R2, t2）。
d. 计算相对位姿：R_rel = R2 * R1.t()，t_rel = t2 - R_rel * t1。
e. 调用cv::stereoRectify：
cv::stereoRectify(
cameraMatrix1, distCoeffs1, cameraMatrix2, distCoeffs2, imageSize, R_rel, t_rel,
R1, R2, P1, P2, Q, flags, alpha, newImageSize, &validPixROI1, &validPixROI2);
注意：这里我们没有畸变系数（distCoeffs1和distCoeffs2设为空）。图像大小是img1.size()。
flags：可以设置为0或CV_CALIB_ZERO_DISPARITY（使两个图像的主点在同一水平线上）。
alpha：我们设置为-1，然后分别得到水平核线影像（校正后的图像）。
返回的R1和R2是旋转矩阵（用于校正），P1和P2是新的投影矩阵，但我们只需要R1和R2来生成校正映射（map）。
f. 使用cv::initUndistortRectifyMap为每个图像生成映射（map1, map2），然后使用cv::remap进行重映射。
g. 对两个图像进行重映射，得到水平核线影像（rectified_img1, rectified_img2）。
h. 竖直核线影像：将原始图像旋转90度，然后重新进行上述校正（注意：旋转90度后，图像大小改变，投影矩阵也需要相应调整？）。或者，我们可以尝试在立体校正时设置不同的参数，但OpenCV不支持竖直核线。因此，我们采用旋转方法：
步骤：
- 将img1和img2分别顺时针旋转90度（这样原来的竖直方向变成了水平方向）。
- 调整投影矩阵：旋转90度相当于对图像应用一个旋转矩阵H_rot（3x3）。在投影矩阵中，我们可以将H_rot乘以P的前3x3部分（即M1和M2），这样新的投影矩阵为P1_rot = H_rot * M1 | H_rot * t1（但注意，平移向量t1部分也要乘以H_rot吗？实际上，投影矩阵的第四列是平移向量，乘以旋转矩阵后，新的投影矩阵为：P1_rot = [H_rot * M1 | H_rot * t1]（因为H_rot是3x3，所以乘以整个3x4矩阵的前三列，第四列是H_rot*t1）。
- 然后，对旋转后的图像和调整后的投影矩阵进行水平核线校正（即上述步骤a-g）。
- 得到校正后的图像后，再逆时针旋转90度，恢复为竖直核线影像（此时核线是竖直的）。
i. 最接近原始图像的核线影像：我们可以尝试不同的alpha值（0.0, 0.5, 1.0）并计算图像变形程度（例如，使用变换矩阵的Frobenius范数减去单位矩阵的范数，或者使用图像特征点匹配的变形程度），然后选择变形最小的。但为了简单，我们可以生成三个不同的结果（alpha=0,0.5,1），然后让用户选择。
16. 由于代码量较大，且需要处理图像旋转和投影矩阵的调整，我们这里只实现水平核线影像和最接近原始图像的核线影像（通过调整alpha）。竖直核线影像的实现较为复杂，我们可能只实现水平核线影像，并说明竖直核线影像可以通过类似方法得到。
17. 考虑到时间，我们简化：只生成水平核线影像，并尝试不同的alpha值（0, 0.5, 1）来生成三个结果，然后选择其中变形最小的作为最接近原始图像的核线影像。如何评价变形？我们可以使用校正变换矩阵H1（或R1）的线性部分（3x3矩阵）与单位矩阵的差异（例如，Frobenius范数）来评价，差异越小，变形越小。
18. 步骤：
a. 对于alpha=0,0.5,1，分别调用stereoRectify，得到R1, R2, P1, P2。
b. 计算每个R1（或H1，这里R1是旋转矩阵，但实际校正变换是H1 = K1_rect * R1 * K1^{-1}）的线性部分（即K1_rect * R1 * K1^{-1}）与单位矩阵的Frobenius范数。选择范数最小的那个alpha对应的结果。
c. 使用该alpha对应的映射进行重映射，得到最接近原始图像的核线影像。
19. 但是，注意：校正变换矩阵H1 = K1_rect * R1 * inv(K1)（其中K1_rect是校正后的内参矩阵，由stereoRectify的P1的前3x3部分给出）。因此，我们需要计算H1，然后计算||H1 - I||_F（Frobenius范数）。但注意，H1通常包含平移和缩放，所以与单位矩阵的差异可能很大。我们也可以只考虑旋转部分，即R1与单位矩阵的差异（因为R1是旋转矩阵，它的范数可以用旋转角度的平方来估计）。
20. 实际上，旋转矩阵R1与单位矩阵的差异可以用旋转角度θ表示：||R1 - I||_F = sqrt(2*(3-2*cosθ-2*cosθ))? 或者更简单地，使用旋转角度θ（通过cv::RQDecomp3x3分解R1得到旋转向量，然后求模）来度量，θ越小越好。
21. 因此，我们可以对每个alpha值，计算R1的旋转角度，选择角度最小的那个。
22. 代码框架：
// 读取图像
Mat img1 = imread("image1.tif", IMREAD_COLOR);
Mat img2 = imread("image2.tif", IMREAD_COLOR);
// 读取投影矩阵文件（假设文件名为projection.txt）
// 文件格式：第一个矩阵3行，每行4个数值；第二个矩阵3行，每行4个数值
// 解析投影矩阵
Mat P1(3,4,CV_64FC1);
Mat P2(3,4,CV_64FC1);
// 分解P1和P2
Mat K1, K2, R1, R2, t1, t2;
decomposeProjectionMatrix(P1, K1, R1, t1); // 注意：OpenCV的decomposeProjectionMatrix函数可以分解投影矩阵
decomposeProjectionMatrix(P2, K2, R2, t2);
// 计算相对位姿
Mat R1_inv = R1.t();
Mat R_rel = R2 * R1_inv;
Mat t_rel = t2 - R_rel * t1;
// 设置相机内参（分解得到的K1和K2是3x3上三角矩阵）
// 注意：decomposeProjectionMatrix分解得到的t是4维的（齐次坐标），我们取前三维
t1 = t1.rowRange(0,3)/t1.at<double>(3);
t2 = t2.rowRange(0,3)/t2.at<double>(3);
// 初始化畸变系数（假设为0）
Mat distCoeffs = Mat::zeros(5,1,CV_64FC1);
// 尝试不同的alpha值
vector<double> alphas = {0.0, 0.5, 1.0};
vector<double> angles; // 存储R1的旋转角度
vector<Mat> rectified_images1; // 存储不同alpha校正后的图像1
vector<Mat> rectified_images2; // 存储不同alpha校正后的图像2
for (double alpha : alphas) {
Mat R1_rect, R2_rect, P1_rect, P2_rect, Q;
Rect validPixROI1, validPixROI2;
stereoRectify(K1, distCoeffs, K2, distCoeffs, img1.size(), R_rel, t_rel,
R1_rect, R2_rect, P1_rect, P2_rect, Q,
CALIB_ZERO_DISPARITY, alpha, img1.size(), &validPixROI1, &validPixROI2);
// 计算R1_rect的旋转角度（将旋转矩阵转换为旋转向量，然后求模）
Vec3d rvec;
Rodrigues(R1_rect, rvec);
double angle = norm(rvec); // 旋转角度（弧度）
angles.push_back(angle);
// 生成映射
Mat map1x, map1y, map2x, map2y;
initUndistortRectifyMap(K1, distCoeffs, R1_rect, P1_rect, img1.size(), CV_32FC1, map1x, map1y);
initUndistortRectifyMap(K2, distCoeffs, R2_rect, P2_rect, img1.size(), CV_32FC1, map2x, map2y);
// 重映射
Mat rect1, rect2;
remap(img1, rect1, map1x, map1y, INTER_LINEAR);
remap(img2, rect2, map2x, map2y, INTER_LINEAR);
rectified_images1.push_back(rect1);
rectified_images2.push_back(rect2);
}
// 选择旋转角度最小的alpha对应的图像
int best_index = min_element(angles.begin(), angles.end()) - angles.begin();
Mat best_rectified_img1 = rectified_images1[best_index];
Mat best_rectified_img2 = rectified_images2[best_index];
// 保存水平核线影像（可以选择任意一个alpha对应的图像，这里我们保存alpha=0.5的作为水平核线影像）
imwrite("rectified_horizontal1.tif", rectified_images1[1]); // alpha=0.5
imwrite("rectified_horizontal2.tif", rectified_images2[1]);
// 保存最接近原始图像的核线影像
imwrite("best_rectified1.tif", best_rectified_img1);
imwrite("best_rectified2.tif", best_rectified_img2);
23. 竖直核线影像：按照上述旋转90度的方法，我们这里不实现，因为代码量很大。如果需要，可以另外实现。
24. 注意：上述代码中，我们假设投影矩阵的分解可以使用OpenCV的`decomposeProjectionMatrix`函数。这个函数可以分解投影矩阵，得到内参矩阵K、旋转矩阵R、平移向量t（以及相机中心）。
25. 由于我们没有具体的图像，所以无法测试。这里提供代码框架。
26. 另外，注意投影矩阵的格式：文件中的矩阵是行优先，每行4个数值。我们读取时按行填充。
27. 完整代码示例（由于没有实际图像，以下代码可能需要调整）：</think>### 核线影像生成代码实现（C++ with OpenCV）
```cpp
#include <opencv2/opencv.hpp>
#include <iostream>
#include <fstream>
#include <vector>
using namespace cv;
using namespace std;
// 读取投影矩阵文件
Mat readProjectionMatrix(const string& filename) {
ifstream file(filename);
if (!file.is_open()) {
cerr << "Error opening: " << filename << endl;
return Mat();
}
string line;
Mat P(3, 4, CV_64F);
int row = 0;
// 跳过"CONTOUR"行
while (getline(file, line)) {
if (line.find("CONTOUR") != string::npos) {
for (int i = 0; i < 3 && getline(file, line); i++) {
istringstream iss(line);
for (int j = 0; j < 4; j++) {
iss >> P.at<double>(i, j);
}
}
break;
}
}
file.close();
return P;
}
// 主函数
int main() {
// 1. 读取影像
Mat img1 = imread("image1.tif", IMREAD_COLOR);
Mat img2 = imread("image2.tif", IMREAD_COLOR);
if (img1.empty() || img2.empty()) {
cerr << "Error loading images!" << endl;
return -1;
}
// 2. 读取投影矩阵
Mat P1 = readProjectionMatrix("proj_matrix1.txt");
Mat P2 = readProjectionMatrix("proj_matrix2.txt");
if (P1.empty() || P2.empty()) {
cerr << "Error loading projection matrices!" << endl;
return -1;
}
// 3. 分解投影矩阵
Mat K1, R1, t1, K2, R2, t2;
decomposeProjectionMatrix(P1, K1, R1, t1);
decomposeProjectionMatrix(P2, K2, R2, t2);
t1 = t1.rowRange(0, 3) / t1.at<double>(3); // 转换为3x1
t2 = t2.rowRange(0, 3) / t2.at<double>(3);
// 4. 计算相对位姿
Mat R_rel = R2 * R1.t();
Mat t_rel = t2 - R_rel * t1;
// 5. 立体校正
Mat R1_rect, R2_rect, P1_rect, P2_rect, Q;
Rect validRoi[2];
stereoRectify(K1, noArray(), K2, noArray(), img1.size(), R_rel, t_rel,
R1_rect, R2_rect, P1_rect, P2_rect, Q,
CALIB_ZERO_DISPARITY, 0, img1.size(), &validRoi[0], &validRoi[1]);
// 6. 生成映射表
Mat map1x, map1y, map2x, map2y;
initUndistortRectifyMap(K1, noArray(), R1_rect, P1_rect, img1.size(), CV_32FC1, map1x, map1y);
initUndistortRectifyMap(K2, noArray(), R2_rect, P2_rect, img2.size(), CV_32FC1, map2x, map2y);
// 7. 水平核线影像
Mat rectified1_h, rectified2_h;
remap(img1, rectified1_h, map1x, map1y, INTER_LINEAR);
remap(img2, rectified2_h, map2x, map2y, INTER_LINEAR);
imwrite("horizontal_epipolar1.tif", rectified1_h);
imwrite("horizontal_epipolar2.tif", rectified2_h);
// 8. 竖直核线影像（通过旋转坐标系）
Mat R_rotate = (Mat_<double>(3,3) << 0, -1, 0, 1, 0, 0, 0, 0, 1); // 90°旋转
Mat P1_rot = R_rotate * P1.colRange(0,3);
hconcat(P1_rot, R_rotate * P1.col(3), P1_rot);
Mat P2_rot = R_rotate * P2.colRange(0,3);
hconcat(P2_rot, R_rotate * P2.col(3), P2_rot);
// 重新分解旋转后的投影矩阵
Mat K1r, R1r, t1r, K2r, R2r, t2r;
decomposeProjectionMatrix(P1_rot, K1r, R1r, t1r);
decomposeProjectionMatrix(P2_rot, K2r, R2r, t2r);
t1r = t1r.rowRange(0,3) / t1r.at<double>(3);
t2r = t2r.rowRange(0,3) / t2r.at<double>(3);
// 竖直方向校正
Mat R1r_rect, R2r_rect, P1r_rect, P2r_rect;
Mat Rr_rel = R2r * R1r.t();
Mat tr_rel = t2r - Rr_rel * t1r;
stereoRectify(K1r, noArray(), K2r, noArray(), img1.size(), Rr_rel, tr_rel,
R1r_rect, R2r_rect, P1r_rect, P2r_rect, Q, CALIB_ZERO_DISPARITY);
// 生成竖直核线影像
Mat map1x_v, map1y_v, map2x_v, map2y_v;
initUndistortRectifyMap(K1r, noArray(), R1r_rect, P1r_rect, img1.size(), CV_32FC1, map1x_v, map1y_v);
initUndistortRectifyMap(K2r, noArray(), R2r_rect, P2r_rect, img2.size(), CV_32FC1, map2x_v, map2y_v);
Mat rectified1_v, rectified2_v;
remap(img1, rectified1_v, map1x_v, map1y_v, INTER_LINEAR);
remap(img2, rectified2_v, map2x_v, map2y_v, INTER_LINEAR);
imwrite("vertical_epipolar1.tif", rectified1_v);
imwrite("vertical_epipolar2.tif", rectified2_v);
// 9. 最接近原始影像的核线影像（最小化几何变形）
vector<double> alphas = {0.0, 0.5, 1.0};
vector<Mat> results;
vector<double> distortions;
for (double alpha : alphas) {
Mat R1a, R2a, P1a, P2a;
stereoRectify(K1, noArray(), K2, noArray(), img1.size(), R_rel, t_rel,
R1a, R2a, P1a, P2a, Q, CALIB_ZERO_DISPARITY, alpha);
// 计算变换矩阵的几何变形
Mat H1 = K1.inv() * R1a * K1;
double distortion = norm(H1 - Mat::eye(3,3,CV_64F), NORM_L2);
distortions.push_back(distortion);
// 生成校正影像
Mat map1x_a, map1y_a;
initUndistortRectifyMap(K1, noArray(), R1a, P1a, img1.size(), CV_32FC1, map1x_a, map1y_a);
Mat rectified;
remap(img1, rectified, map1x_a, map1y_a, INTER_LINEAR);
results.push_back(rectified);
}
// 选择最小变形的结果
int best_idx = min_element(distortions.begin(), distortions.end()) - distortions.begin();
imwrite("minimal_distortion.tif", results[best_idx]);
cout << "核线影像生成完成！" << endl;
return 0;
}
```
### 关键步骤说明：
1. **投影矩阵解析**：读取包含`CONTOUR`标记的文本文件，提取3×4投影矩阵
2. **矩阵分解**：使用`decomposeProjectionMatrix`分解内参(K)、旋转(R)和平移(t)
3. **立体校正**：
- 水平核线：直接计算校正映射，生成水平对齐影像
- 竖直核线：通过坐标系旋转（90°）实现核线方向转换
4. **最小变形处理**：通过调整`stereoRectify`的alpha参数(0.0/0.5/1.0)，选择几何变形最小的结果
### 输出文件：
- `horizontal_epipolar[1|2].tif`：水平核线影像
- `vertical_epipolar[1|2].tif`：竖直核线影像
- `minimal_distortion.tif`：最接近原始影像的核线结果
本文介绍了一种求解不同最大子矩阵个数的高效算法，通过悬线法结合哈希技术实现矩阵的快速去重，确保计算准确性和效率。
扫一扫
3
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
