#include <opencv2/opencv.hpp>
#include <iostream>
#include <vector>
using namespace cv;
using namespace std;
class DigitalDesign {
public:
// 计算相机中心
Mat getCameraCenter(const Mat& P) {
Mat C;
Mat R = P.colRange(0, 3);
Mat T = P.col(3);
Mat inv_R;
invert(R, inv_R);
C = -1 * inv_R * T;
return C;
}
// 计算相机焦距
double getCameraf(const Mat& K) {
double f;
f = (abs(K.at<double>(0, 0)) + abs(K.at<double>(1, 1))) / 2;
return f;
}
// 计算单应矩阵实现核线校正
Mat computeHomography(const Mat& K, const Mat& R, const Mat& rop, const Mat& kop) {
// 单应矩阵 H = K * R * rop^(-1) * kop^(-1)
return K * R * rop.t() * kop.inv();
}
// 计算角点,确定影像范围
void getCornerpoints(const Mat& R, const Mat& K, const Mat& K_temp, const Mat& rop, const Mat& imsize,
vector<double>& w_min, vector<double>& w_max, vector<double>& h_min, vector<double>& h_max) {
Mat corner_Pts(3, 4, CV_64FC1); // 图像四个角点
corner_Pts.at<double>(0, 0) = 0;
corner_Pts.at<double>(1, 0) = 0;
corner_Pts.at<double>(2, 0) = 1;
corner_Pts.at<double>(0, 1) = imsize.at<double>(0, 1) - 1;
corner_Pts.at<double>(1, 1) = 0;
corner_Pts.at<double>(2, 1) = 1;
corner_Pts.at<double>(0, 2) = 0;
corner_Pts.at<double>(1, 2) = imsize.at<double>(0, 0) - 1;
corner_Pts.at<double>(2, 2) = 1;
corner_Pts.at<double>(0, 3) = imsize.at<double>(0, 1) - 1;
corner_Pts.at<double>(1, 3) = imsize.at<double>(0, 0) - 1;
corner_Pts.at<double>(2, 3) = 1;
Mat inv_K;
invert(K, inv_K);
Mat epi_corner_Pts = K_temp * rop * R.t() * inv_K * corner_Pts;
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 3; ++j) {
epi_corner_Pts.at<double>(j, i) /= epi_corner_Pts.at<double>(2, i); // 齐次化
}
}
double min_w, max_w;
minMaxLoc(epi_corner_Pts.row(0), &min_w, &max_w);
double min_h, max_h;
minMaxLoc(epi_corner_Pts.row(1), &min_h, &max_h);
w_min.push_back(min_w);
w_max.push_back(max_w);
h_min.push_back(min_h);
h_max.push_back(max_h);
}
// 生成核线影像(使用单应矩阵)
Mat generateEpipolarLineImagesWithHomography(const Mat& image, const Mat& H, const Size& size) {
Mat epi_image;
// 使用单应矩阵进行透视变换
warpPerspective(image, epi_image, H.inv(), size);
return epi_image;
}
};
int main() {
Mat image1 = imread("015AR027.tif", IMREAD_COLOR);
Mat image2 = imread("015AR028.tif", IMREAD_COLOR);
if (image1.empty() || image2.empty()) {
cout << "无法读取图像!请检查文件路径是否正确。" << endl;
return -1;
}
// 创建投影矩阵P
Mat P1 = (Mat_<double>(3, 4) << -12179.8753230552230, 7050.9952897638232, -2467.7527210258945, -17321440031.6463700,
4723.7392252925192, 5779.0243267188052, -11880.0042304789190, -22040783539.28058600,
-0.4907881086633, -0.5142398041434, -0.7033380810316, 2000536.13817580420);
Mat P2 = (Mat_<double>(3, 4) <<
-12158.6697497674500, 7092.3578385200199, -2453.7545028069599, -17474815958.9305270000000,
4621.9021026258170, 5524.7193480639962, -12039.9681076173220, -21127026467.7189250000000,
-0.5027584877552, -0.5235467256044, -0.6878464429645, 2038373.8802346450000);
// 定义内参矩阵K、旋转矩阵R和位移矩阵T
Mat K1, R1, t1; Mat K2, R2, t2;
// 分解投影矩阵P->K,R,T
decomposeProjectionMatrix(P1, K1, R1, t1);
// R1化为正
double det_R1 = determinant(R1);
if (det_R1 < 0) {
R1 = -R1;
}
// 去除齐次坐标
double tt1 = t1.at<double>(3, 0);
Mat T1 = t1 / tt1;
T1 = T1.rowRange(0, 3);
decomposeProjectionMatrix(P2, K2, R2, t2);
// R2化为正
double det_R2 = determinant(R2);
if (det_R2 < 0) {
R2 = -R2;
}
// 去除齐次坐标
double tt2 = t2.at<double>(3, 0);
Mat T2 = t2 / tt2;
T2 = T2.rowRange(0, 3);
// 创建DigitalDesign类的实例
DigitalDesign dd;
// 相机中心坐标
Mat C1 = dd.getCameraCenter(P1);
Mat C2 = dd.getCameraCenter(P2);
// 相机的焦距
double f1 = dd.getCameraf(K1);
double f2 = dd.getCameraf(K2);
// 计算基线
Mat B = C2 - C1;
// 设置核线影像像平面(水平方向)
Mat temp = (Mat_<double>(3, 1) << 0, 0, 1); temp = temp.t();
Mat r1, r2, r3;
r1 = B / norm(B); r1 = r1.t();
r2 = temp.cross(B.t()) / norm(temp.cross(B.t()));
r3 = r1.cross(r2) / norm(r1.cross(r2));
// 核线影像旋转矩阵 rop
Mat rop(3, 3, CV_64FC1);
r1.copyTo(rop.row(0));
r2.copyTo(rop.row(1));
r3.copyTo(rop.row(2));
// 计算核线影像焦距
double kop_f;
Mat F1 = abs(rop.row(2) * R1.row(2).t() * f1);
Mat F2 = abs(rop.row(2) * R2.row(2).t() * f2);
kop_f = (F1.at<double>(0, 0) > F2.at<double>(0, 0)) ? F2.at<double>(0, 0) : F1.at<double>(0, 0);
// 像主点为0时的内矩阵
Mat K_temp = (Mat_<double>(3, 3) << -kop_f, 0, 0, 0, -kop_f, 0, 0, 0, 1);
// 图像尺寸信息
Mat imsize1(1, 3, CV_64FC1);
imsize1.at<double>(0, 0) = image1.rows; // 高度
imsize1.at<double>(0, 1) = image1.cols; // 宽度
imsize1.at<double>(0, 2) = image1.channels();
Mat imsize2(1, 3, CV_64FC1);
imsize2.at<double>(0, 0) = image2.rows;
imsize2.at<double>(0, 1) = image2.cols;
imsize2.at<double>(0, 2) = image2.channels();
// 计算角点,确定影像范围
vector<double> w_min, w_max, h_min, h_max;
dd.getCornerpoints(R1, K1, K_temp, rop, imsize1, w_min, w_max, h_min, h_max);
dd.getCornerpoints(R2, K2, K_temp, rop, imsize2, w_min, w_max, h_min, h_max);
int epi_h[2] = { round(h_max[0] - h_min[0]), round(h_max[1] - h_min[1]) };
int epi_w[2] = { round(w_max[0] - w_min[0]), round(w_max[1] - w_min[1]) };
// 核线影像的像主点
Mat PP1 = (Mat_<double>(1, 2) << -w_min[0], h_max[0]);
Mat PP2 = (Mat_<double>(1, 2) << -w_min[1], h_max[1]);
// 计算核线影像的内参矩阵 kop
Mat kop1 = (Mat_<double>(3, 3) << -kop_f, 0, PP1.at<double>(0, 0),
0, kop_f, PP1.at<double>(0, 1),
0, 0, 1);
Mat kop2 = (Mat_<double>(3, 3) << -kop_f, 0, PP2.at<double>(0, 0),
0, kop_f, PP2.at<double>(0, 1),
0, 0, 1);
// 计算单应矩阵
Mat H1 = dd.computeHomography(K1, R1, rop, kop1);
Mat H2 = dd.computeHomography(K2, R2, rop, kop2);
// 确定输出图像尺寸(取两者最大值以包含所有内容)
Size epi_size(max(epi_w[0], epi_w[1]), max(epi_h[0], epi_h[1]));
// 生成核线影像(使用单应矩阵)
Mat epi_image1 = dd.generateEpipolarLineImagesWithHomography(image1, H1, epi_size);
Mat epi_image2 = dd.generateEpipolarLineImagesWithHomography(image2, H2, epi_size);
// 保存结果(请替换为你的保存路径)
imwrite("EpiImage_Left_Homography.jpg", epi_image1);
imwrite("EpiImage_Right_Homography.jpg", epi_image2);
cout << "核线影像已保存。" << endl;
// 显示图像
namedWindow("原始左图", WINDOW_NORMAL);
namedWindow("原始右图", WINDOW_NORMAL);
namedWindow("核线左图", WINDOW_NORMAL);
namedWindow("核线右图", WINDOW_NORMAL);
imshow("原始左图", image1);
imshow("原始右图", image2);
imshow("核线左图", epi_image1);
imshow("核线右图", epi_image2);
// 显示匹配效果(在核线影像上绘制水平线)
for (int y = 0; y < epi_size.height; y += 50) {
line(epi_image1, Point(0, y), Point(epi_size.width, y), Scalar(0, 255, 0), 1);
line(epi_image2, Point(0, y), Point(epi_size.width, y), Scalar(0, 255, 0), 1);
}
namedWindow("核线左图(带水平线)", WINDOW_NORMAL);
namedWindow("核线右图(带水平线)", WINDOW_NORMAL);
imshow("核线左图(带水平线)", epi_image1);
imshow("核线右图(带水平线)", epi_image2);
waitKey(0);
return 0;
}把rop中的temp改为基线和(0,0,1)和基线的叉乘实现竖核线影像生成以及temp变为r1实现最接近原始图像的核线影像