跟踪预测后,已经能建立多帧之间标牌的一致性关系,而每一帧的外参通过前面的定位环节已经确定了,因此只要建立同名点关系,就可以用三角量测把坐标算出来了。
有一种办法是采用特征匹配的方法,通过特征点匹配建立同名点关系,而标牌还多一个单应约束,建立最小二乘问题,问题即可解。
但本文并不采用这种方法,特征匹配强依赖于特征点的个数,对照片模糊或特征比较少的要素就无能为力了,本文介绍一种更加通用的方法,直接使用目标检测框来做:
基于这样的判断:标牌角点投影到各帧图像中与目标检测的角点重投影误差时的位置就是标牌的位置,基于此建立非线性优化算法,问题得解。
算法步骤如下:
1.每一帧两两计算相机中心与标牌中心构成的射线的夹角,选出其中最大夹角,这一步是因为三角量算的精度跟基线长度是呈正比关系,因此最大夹角就能保证基线最长.
2、将第一步选出的两帧,将标牌中心当成同名点三角化得到三维坐标:
如何三角化同名点,我们知道三维点(X,Y,Z)投影到图像平面上是通过:
其中,是一个3*4的矩阵,设为P,则有:
两边同时叉乘u
展开后有:
继续代入得:
上式三个方程只能提供两个约束,因为
每一帧都可以提供两个约束,X总共有三个未知数,二帧就可以求解最小二乘算法了,代码如下:
bool SignVec::triangulateView(const SignVecMeta& to_vec1, const SignVecMeta& to_vec2, const Camera& camera, Eigen::Vector3d& vec3d)
{
Eigen::Matrix3d cRb = camera.mount.rotation.matrix();
Eigen::Vector3d cTb = camera.mount.translation;
Eigen::Matrix3d bRw1 = to_vec1.pose.orientation.toRotationMatrix();
Eigen::Vector3d bTw1 = to_vec1.pose.position;
Eigen::Matrix3d R1 = bRw1 * cRb;
Eigen::Vector3d C1 = bTw1 + bRw1 * cTb;
Eigen::Matrix<double, 3, 4> P1 = Eigen::Matrix<double, 3, 4>::Zero();
P1.block(0, 0, 3, 3) = R1.transpose();
P1.col(3) = - R1.transpose() * C1;
P1 = camera.K * P1;
Eigen::Matrix3d bRw2 = to_vec2.pose.orientation.toRotationMatrix();
Eigen::Vector3d bTw2 = to_vec2.pose.position;
Eigen::Matrix3d R2 = bRw2 * cRb;
Eigen::Vector3d C2 = bTw2 + bRw2 * cTb;
Eigen::Matrix<double, 3, 4> P2 = Eigen::Matrix<double, 3, 4>::Zero();
P2.block(0, 0, 3, 3) = R2.transpose();
P2.col(3) = - R2.transpose() * C2;
P2 = camera.K * P2;
const Eigen::Matrix<double, 1, 4>& P11 = P1.block(0, 0, 1, 4);
const Eigen::Matrix<double, 1, 4>& P12 = P1.block(1, 0, 1, 4);
const Eigen::Matrix<double, 1, 4>& P13 = P1.block(2, 0, 1, 4);
const Eigen::Matrix<double, 1, 4>& P21 = P2.block(0, 0, 1, 4);
const Eigen::Matrix<double, 1, 4>& P22 = P2.block(1, 0, 1, 4);
const Eigen::Matrix<double, 1, 4>& P23 = P2.block(2, 0, 1, 4);
Eigen::Vector2d cent2d1 = to_vec1.dlSign.bBox.GetCenter();
Eigen::Vector2d cent2d2 = to_vec2.dlSign.bBox.GetCenter();
const double &u1 = cent2d1.x();
const double &v1 = cent2d1.y();
const double &u2 = cent2d2.x();
const double &v2 = cent2d2.y();
Eigen::Matrix4d H;
H.block(0, 0, 1, 4) = v1 * P13 - P12;
H.block(1, 0, 1, 4) = P11 - u1 * P13;
H.block(2, 0, 1, 4) = v2 * P23 - P22;
H.block(3, 0, 1, 4) = P21 - u2 * P23;
// SVD
Eigen::JacobiSVD<Eigen::MatrixXd> svd(H, Eigen::ComputeFullU | Eigen::ComputeFullV);
Eigen::Matrix4d V = svd.matrixV();
Eigen::Vector4d X = V.block(0, 3, 4, 1);
X = X / X(3, 0);
vec3d = X.block(0, 0, 3, 1);
return true;
}
3.将第2步得到的三维坐标反投到第一帧上,作为深度初始值,这将代入第4步进行非线性优化:
double SignVec::calculateDepth(const Pose& pose, const Camera& camera, Eigen::Vector3d& vec3d) const
{
Eigen::Matrix3d cRb = camera.mount.rotation.matrix();
Eigen::Vector3d cTb = camera.mount.translation;
Eigen::Matrix3d bRw0 = pose.orientation.toRotationMatrix();
Eigen::Vector3d bTw0 = pose.position;
Eigen::Matrix3d R0 = bRw0 * cRb;
Eigen::Vector3d C0 = bTw0 + bRw0 * cTb;
Eigen::Vector3d camera3D = R0.transpose() * (vec3d - C0);
return camera3D.z();
}
4.将第二步得到深度将第一帧的标牌角点投影到三维坐标系中。将逆深度作为优化量,将标牌角点到每一帧图像的反投影误差作为残差,建立非线性优化问题。(为了使其更符合正态分布,这里采用的是逆深度,即深度的导数)
构建残差:
因此:
代码如下:
void SignVec::refineDepth(const std::vector<SignVecMeta> &to_vecs, const Camera &camera, double &depth) const
{
Eigen::Matrix3d cRb = camera.mount.rotation.matrix();
Eigen::Vector3d cTb = camera.mount.translation;
int N = to_vecs.size();
std::vector<Eigen::Vector2d> q1(4), q2(4);
q1[0] = Eigen::Vector2d(to_vecs[0].dlSign.bBox.GetMinx(), to_vecs[0].dlSign.bBox.GetMiny());
q1[1] = Eigen::Vector2d(to_vecs[0].dlSign.bBox.GetMinx() + to_vecs[0].dlSign.bBox.Width(), to_vecs[0].dlSign.bBox.GetMiny());
q1[2] = Eigen::Vector2d(to_vecs[0].dlSign.bBox.GetMinx() + to_vecs[0].dlSign.bBox.Width(), to_vecs[0].dlSign.bBox.GetMiny() + to_vecs[0].dlSign.bBox.Height());
q1[3] = Eigen::Vector2d(to_vecs[0].dlSign.bBox.GetMinx(), to_vecs[0].dlSign.bBox.GetMiny() + to_vecs[0].dlSign.bBox.Height());
Eigen::Matrix3d bRw0 = to_vecs[0].pose.orientation.toRotationMatrix();
Eigen::Vector3d bTw0 = to_vecs[0].pose.position;
Eigen::Matrix3d R0 = bRw0 * cRb;
Eigen::Vector3d C0 = bTw0 + bRw0 * cTb;
double inverse_depth = 1.0 / depth;
std::vector<Eigen::Vector3d> Points3d(4), Ray(4);
for (int i = 0; i < 4; i++)
{
//let ray length aligned with body is as same as depth
Ray[i] = cRb * camera.K.inverse() * Eigen::Vector3d(q1[i].x(), q1[i].y(), 1.0);
Ray[i] = Ray[i] / Ray[i].x();
Ray[i] = bRw0 * Ray[i];
Ray[i].normalize();
Points3d[i] = C0 + depth * Ray[i];
}
double stopThresholdLM = 0.0;
double currentEro = 0.0;
double currentLambda = 0.0;
double ni = 2.0;
//Generate init lamda
{
Eigen::MatrixXd H = Eigen::MatrixXd::Zero(1, 1);
for(int i = 1; i < N; i++)
{
q2[0] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx(), to_vecs[i].dlSign.bBox.GetMiny());
q2[1] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx() + to_vecs[i].dlSign.bBox.Width(), to_vecs[i].dlSign.bBox.GetMiny());
q2[2] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx() + to_vecs[i].dlSign.bBox.Width(), to_vecs[i].dlSign.bBox.GetMiny() + to_vecs[i].dlSign.bBox.Height());
q2[3] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx(), to_vecs[i].dlSign.bBox.GetMiny() + to_vecs[i].dlSign.bBox.Height());
Eigen::Matrix3d bRw = to_vecs[i].pose.orientation.toRotationMatrix();
Eigen::Vector3d bTw = to_vecs[i].pose.position;
Eigen::Matrix3d R = bRw * cRb;
Eigen::Vector3d C = bTw + bRw * cTb;
for (int j = 0; j < 4; j++)
{
Eigen::Vector3d camera3D = R.transpose() * (Points3d[j] - C);
if (camera3D[2] <= 0.0)
continue;
Eigen::Vector3d uvw = camera.K * camera3D;
Eigen::Vector2d uv(uvw[0] / uvw[2], uvw[1] / uvw[2]);
Eigen::Vector2d residual = uv - q2[j];
currentEro += residual.transpose() * residual;
Eigen::Matrix<double, 2, 3> J_F;
J_F << camera.fx() / camera3D.z(), 0.0, -camera.fx() * camera3D.x() / (camera3D.z() * camera3D.z()),
0.0, camera.fy() / camera3D.z(), -camera.fy() * camera3D.y() / (camera3D.z() * camera3D.z());
Eigen::Matrix<double, 2, 1> J_aux = J_F * R.transpose() * (- 1.0 / (inverse_depth * inverse_depth)) * Ray[j];
H += J_aux.transpose() * J_aux;
}
}
stopThresholdLM = 1e-6 * currentEro;
double maxDiagonal = std::max(std::fabs(H(0, 0)), maxDiagonal);
double tau = 1e-5;
currentLambda = tau * maxDiagonal;
}
auto robustWeightCauchy = [=](double norm_res) -> double {
return std::sqrt(1.0 / (1.0 + norm_res * norm_res));
};
int iter = 0;
bool stop = false;
while (!stop && iter < 30)
{
bool oneStepSuccess = false;
int false_cnt = 0;
while (!oneStepSuccess)
{
Eigen::MatrixXd H = Eigen::MatrixXd::Zero(1, 1);
Eigen::MatrixXd g = Eigen::MatrixXd::Zero(1, 1);
for (int i = 1; i < N; i++)
{
Points3d[0] = C0 + 1.0 / inverse_depth * Ray[0];
Points3d[1] = C0 + 1.0 / inverse_depth * Ray[1];
Points3d[2] = C0 + 1.0 / inverse_depth * Ray[2];
Points3d[3] = C0 + 1.0 / inverse_depth * Ray[3];
q2[0] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx(), to_vecs[i].dlSign.bBox.GetMiny());
q2[1] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx() + to_vecs[i].dlSign.bBox.Width(), to_vecs[i].dlSign.bBox.GetMiny());
q2[2] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx() + to_vecs[i].dlSign.bBox.Width(), to_vecs[i].dlSign.bBox.GetMiny() + to_vecs[i].dlSign.bBox.Height());
q2[3] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx(), to_vecs[i].dlSign.bBox.GetMiny() + to_vecs[i].dlSign.bBox.Height());
Eigen::Matrix3d bRw = to_vecs[i].pose.orientation.toRotationMatrix();
Eigen::Vector3d bTw = to_vecs[i].pose.position;
Eigen::Matrix3d R = bRw * cRb;
Eigen::Vector3d C = bTw + bRw * cTb;
for (int j = 0; j < 4; j++)
{
Eigen::Vector3d camera3D = R.transpose() * (Points3d[j] - C);
if (camera3D[2] <= 0.0)
continue;
Eigen::Vector3d uvw = camera.K * camera3D;
Eigen::Vector2d uv(uvw[0] / uvw[2], uvw[1] / uvw[2]);
Eigen::Vector2d residual = uv - q2[j];
double err_i_norm = residual.norm();
double weight = robustWeightCauchy(err_i_norm);
Eigen::Matrix<double, 2, 3> J_F;
J_F << camera.fx() / camera3D.z(), 0.0, -camera.fx() * camera3D.x() / (camera3D.z() * camera3D.z()),
0.0, camera.fy() / camera3D.z(), -camera.fy() * camera3D.y() / (camera3D.z() * camera3D.z());
Eigen::Matrix<double, 2, 1> J_aux = J_F * R.transpose() * (- 1.0 / (inverse_depth * inverse_depth)) * Ray[j];
H += J_aux.transpose() * J_aux * weight;
g -= J_aux.transpose() * residual * weight;
}
}
H(0, 0) += currentLambda;
Eigen::Matrix<double, 1, 1> delta = H.colPivHouseholderQr().solve(g);
if (delta.squaredNorm() < EPS || false_cnt > 10)
{
stop = true;
break;
}
inverse_depth += delta[0];
double scale = 0;
scale = delta.transpose() * (currentLambda * delta + g);
scale += 1e-3;
double temp_err = 0.0;
for (int i = 1; i < N; i++)
{
Points3d[0] = C0 + 1.0 / inverse_depth * Ray[0];
Points3d[1] = C0 + 1.0 / inverse_depth * Ray[1];
Points3d[2] = C0 + 1.0 / inverse_depth * Ray[2];
Points3d[3] = C0 + 1.0 / inverse_depth * Ray[3];
q2[0] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx(), to_vecs[i].dlSign.bBox.GetMiny());
q2[1] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx() + to_vecs[i].dlSign.bBox.Width(), to_vecs[i].dlSign.bBox.GetMiny());
q2[2] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx() + to_vecs[i].dlSign.bBox.Width(), to_vecs[i].dlSign.bBox.GetMiny() + to_vecs[i].dlSign.bBox.Height());
q2[3] = Eigen::Vector2d(to_vecs[i].dlSign.bBox.GetMinx(), to_vecs[i].dlSign.bBox.GetMiny() + to_vecs[i].dlSign.bBox.Height());
Eigen::Matrix3d bRw = to_vecs[i].pose.orientation.toRotationMatrix();
Eigen::Vector3d bTw = to_vecs[i].pose.position;
Eigen::Matrix3d R = bRw * cRb;
Eigen::Vector3d C = bTw + bRw * cTb;
for (int j = 0; j < 4; j++)
{
Eigen::Vector3d camera3D = R.transpose() * (Points3d[j] - C);
if (camera3D[2] <= 0.0)
continue;
Eigen::Vector3d uvw = camera.K * camera3D;
Eigen::Vector2d uv(uvw[0] / uvw[2], uvw[1] / uvw[2]);
Eigen::Vector2d residual = uv - q2[j];
temp_err += residual.transpose() * residual;
}
}
double rho = (currentEro - temp_err) / scale;
if (rho > 0 && std::isfinite(temp_err)) // last step was good
{
double alpha = 1.0 - pow((2 * rho - 1), 3);
alpha = std::min(alpha, 2. / 3.);
double scaleFactor = (std::max)(1. / 3., alpha);
currentLambda *= scaleFactor;
ni = 2;
currentEro = temp_err;
oneStepSuccess = true;
}
else
{
currentLambda *= ni;
ni *= 2;
oneStepSuccess = false;
}
if(!oneStepSuccess)
{
depth -= delta[0];
false_cnt++;
}
else
{
false_cnt = 0;
}
}
iter++;
if (sqrt(currentEro) <= stopThresholdLM)
stop = true;
}
depth = 1.0 / inverse_depth;
}
最后的重建结果如图: