MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
homography.cc
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#include <set>
11#include <stdexcept>
12#include <vector>
13
14#include "math/matrix_tools.h"
15#include "math/matrix_svd.h"
16#include "sfm/homography.h"
17
19
20bool
22{
23 if (points.size() < 4)
24 throw std::invalid_argument("At least 4 matches required");
25
26 /* Create 2Nx9 matrix A. Each correspondence creates two rows in A. */
27 std::vector<double> A(2 * points.size() * 9);
28 for (std::size_t i = 0; i < points.size(); ++i)
29 {
30 std::size_t const row1 = 9 * i;
31 std::size_t const row2 = 9 * (i + points.size());
32 Correspondence2D2D const& match = points[i];
33 A[row1 + 0] = 0.0;
34 A[row1 + 1] = 0.0;
35 A[row1 + 2] = 0.0;
36 A[row1 + 3] = match.p1[0];
37 A[row1 + 4] = match.p1[1];
38 A[row1 + 5] = 1.0;
39 A[row1 + 6] = -match.p1[0] * match.p2[1];
40 A[row1 + 7] = -match.p1[1] * match.p2[1];
41 A[row1 + 8] = -match.p2[1];
42 A[row2 + 0] = -match.p1[0];
43 A[row2 + 1] = -match.p1[1];
44 A[row2 + 2] = -1.0;
45 A[row2 + 3] = 0.0;
46 A[row2 + 4] = 0.0;
47 A[row2 + 5] = 0.0;
48 A[row2 + 6] = match.p1[0] * match.p2[0];
49 A[row2 + 7] = match.p1[1] * match.p2[0];
50 A[row2 + 8] = match.p2[0];
51 }
52
53 /* Compute homography matrix using SVD. */
55 math::matrix_svd<double>(&A[0], 2 * points.size(),
56 9, nullptr, nullptr, V.begin());
57
58 /* Only consider the last column of V as the solution. */
59 for (int i = 0; i < 9; ++i)
60 (*result)[i] = V[i * 9 + 8];
61
62 return true;
63}
64
65double
67 Correspondence2D2D const& match)
68{
69 /*
70 * Computes the symmetric transfer error for a given match and homography
71 * matrix. The error is computed as [Sect 4.2.2, Hartley, Zisserman]:
72 *
73 * e = d(x, (H^-1)x')^2 + d(x', Hx)^2
74 */
75 math::Vec3d p1(match.p1[0], match.p1[1], 1.0);
76 math::Vec3d p2(match.p2[0], match.p2[1], 1.0);
77
78 math::Matrix3d invH = math::matrix_inverse(homography);
79 math::Vec3d result = invH * p2;
80 result /= result[2];
81 double error = (p1 - result).square_norm();
82
83 result = homography * p1;
84 result /= result[2];
85 error += (result - p2).square_norm();
86
87 return 0.5 * error;
88}
89
T * begin(void)
Definition matrix.h:506
Matrix< T, N, N > matrix_inverse(Matrix< T, N, N > const &mat)
Calculates the inverse of the given matrix.
bool homography_dlt(Correspondences2D2D const &points, HomographyMatrix *result)
Direct linear transformation algorithm to compute the homography matrix from image correspondences.
Definition homography.cc:21
std::vector< Correspondence2D2D > Correspondences2D2D
double symmetric_transfer_error(HomographyMatrix const &homography, Correspondence2D2D const &match)
Computes the symmetric transfer error for an image correspondence given the homography matrix between...
Definition homography.cc:66
#define SFM_NAMESPACE_END
Definition defines.h:14
#define SFM_NAMESPACE_BEGIN
Definition defines.h:13
Two image coordinates which correspond to each other in terms of observing the same point in the scen...