-
Notifications
You must be signed in to change notification settings - Fork 71
/
FMatrixAffineCompute.cpp
65 lines (43 loc) · 1.38 KB
/
FMatrixAffineCompute.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include <vnl/algo/vnl_svd.h>
template <class T>
T sum(vector<T> const& vector) {
T sum = vector[0];
for(int i=1; i<vector.size(); i++)
sum += vector[i];
return sum;
}
Vector computeAffineFMatrix(vector<HomPoint> const& points1,
vector<HomPoint> const& points2)
{
assert(points1.size() == points2.size());
int n = points1.size();
Vector centroid(4, 0.0);
for(int i=0; i<n; i++) {
centroid[0] += points1[i].x();
centroid[1] += points1[i].y();
centroid[2] += points2[i].x();
centroid[3] += points2[i].y();
}
centroid /= n;
Matrix matrix(n, 4);
for(int i=0; i<n; i++) {
matrix(i, 0) = points1[i].x() - centroid[0];
matrix(i, 1) = points1[i].y() - centroid[1];
matrix(i, 2) = points2[i].x() - centroid[2];
matrix(i, 3) = points2[i].y() - centroid[3];
}
// cout << "fmatrixcompute: " << endl << matrix << endl << endl;
vnl_svd<double> svd(matrix);
// cout << "U = " << endl << svd.U() << endl << endl;
// cout << "W = " << endl << svd.W() << endl << endl;
// cout << "V = " << endl << svd.V() << endl << endl;
Vector v = svd.V().get_column(3);
Vector result(5);
result.update(v); // a,b,c,d
result[4] = inner_product(v, centroid);
if (result[4] < 0)
result = -result;
// cout << "mat: " << v << endl;
// cout << "test: " << matrix * v << endl;
return result.normalize();
}