从两幅图像的匹配点计算焦距f

从网络中找到两种计算焦距的方法
1。两幅图像的单应矩阵H计算焦距f（取自《全景图拼接》opencv代码）

//从两幅图像的单应矩阵H计算焦距f。
void focalsFromHomography(const Mat& H, double &f0, double &f1, bool &f0_ok, bool &f1_ok)
{
    CV_Assert(H.type() == CV_64F && H.size() == Size(3, 3));

    const double* h = H.ptr<double>();//h 是一个向量

    double d1, d2; // 分母
    double v1, v2; // 焦距平方值候选

    f1_ok = true;
    d1 = h[6] * h[7];
    d2 = (h[7] - h[6]) * (h[7] + h[6]);
    v1 = -(h[0] * h[1] + h[3] * h[4]) / d1;
    v2 = (h[0] * h[0] + h[3] * h[3] - h[1] * h[1] - h[4] * h[4]) / d2;
    if (v1 < v2) std::swap(v1, v2);
    if (v1 > 0 && v2 > 0) f1 = std::sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
    else if (v1 > 0) f1 = std::sqrt(v1);
    else f1_ok = false;

    f0_ok = true;
    d1 = h[0] * h[3] + h[1] * h[4];
    d2 = h[0] * h[0] + h[1] * h[1] - h[3] * h[3] - h[4] * h[4];
    v1 = -h[2] * h[5] / d1;
    v2 = (h[5] * h[5] - h[2] * h[2]) / d2;
    if (v1 < v2) std::swap(v1, v2);
    if (v1 > 0 && v2 > 0) f0 = std::sqrt(std::abs(d1) > std::abs(d2) ? v1 : v2);
    else if (v1 > 0) f0 = std::sqrt(v1);
    else f0_ok = false;
}

用下面的测试

	//读入两个图
	Mat img1 = imread(name1);//
	Mat img2 = imread(name2);//

	//匹配
	vector<KeyPoint> kp1,kp2;
	vector<DMatch> matches;
	meMatch2(img1, img2,kp1,kp2,matches);


	//生成匹配点
	vector<Point2f> points1, points2;
	get_match_points (kp1, kp2, matches, points1, points2);

	//4对相对应的像点计算图像的单应 Homography
	Mat H = findHomography(points1, points2, 0);

	double f0; double f1; bool f0_ok; bool f1_ok;
	//从两幅图像的单应矩阵H计算焦距f。
	focalsFromHomography(H, f0, f1, f0_ok, f1_ok);
	if (f0_ok)
			std::cout <<"图1焦距："<< f0 << '\n';
	if (f1_ok)
			std::cout <<"图2焦距："<< f1 << '\n';

结果：

图1焦距：81.9963
图2焦距：992.465
只有图2计算得很算准确。

-=-=-=-=-=-=-=-=-=-----------------分隔线------------------------=-=-=-=-=-=-=-=-=

2。两图的匹配点计算焦距f(取自一个《极线校正工具包》matlab)

function [f] = compIntrinsic (ml,mr,width,height)
%计算一个焦距
%
%
% ml,mr 是图像中的对应点 
%
% width, height 是图像的尺寸(用于取得焦距和图像中心)
% 



options = optimset('Jacobian','off');

a0 = [0 0 0 0 0, 0];
[af,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@(x) costRectif(x,width,height,ml,mr), a0,...
    [-pi/2, -pi/2, -pi/2, -pi/2, -pi/2, -1.5], [pi/2, pi/2, pi/2, pi/2, pi/2, 1.5],options);

counter = 0;
while af(6)>1 || af(6)<-1
    counter = counter+1;

    disp('焦点超出范围；正在重新尝试(随机)...');

    a0 = [0 0 0 0 0, 2*rand-1.0];

    [af,resnorm] = lsqnonlin(@(x) costRectif(x,width,height,ml,mr), a0,...
        [-pi/2, -pi/2, -pi/2, -pi/2, -pi/2, -1.5], [pi/2, pi/2, pi/2, pi/2, pi/2, 1.5],options);

    if counter > 20
        warning('找不到范围内的焦点')
        break
    end

end

f =3^af(6)*(width+height);
% disp(['焦距：f = ' num2str(f)]);

测试

f = compIntrinsic (ml,mr,width,height);%ml,mr 一个2xn矩阵的匹配点对.width,height 是图像的长，宽
disp(['焦距：f = ' num2str(f)]);

结果
焦距：f = 892.7031
基本准确。