← Back to branch summary

~siretart/ubuntu/utopic/blender/libav10

« back to all changes in this revision

Viewing changes to extern/libmv/libmv/simple_pipeline/initialize_reconstruction.cc

  • Committer: Package Import Robot
  • Author(s): Matteo F. Vescovi
  • Date: 2012-07-23 08:54:18 UTC
  • mfrom: (14.2.16 sid)
  • mto: (14.2.19 sid)
  • mto: This revision was merged to the branch mainline in revision 42.
  • Revision ID: package-import@ubuntu.com-20120723085418-9foz30v6afaf5ffs
Tags: 2.63a-2
http://bugs.debian.org/658075
* debian/: Cycles support added (Closes: #658075)
  For now, this top feature has been enabled only
  on [any-amd64 any-i386] architectures because
  of OpenImageIO failing on all others
* debian/: scripts installation path changed
  from /usr/lib to /usr/share:
  + debian/patches/: patchset re-worked for path changing
  + debian/control: "Breaks" field added on yafaray-exporter

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/libmv/libmv/simple_pipeline/initialize_reconstruction.cc
 
1
// Copyright (c) 2011 libmv authors.
 
2
// 
 
3
// Permission is hereby granted, free of charge, to any person obtaining a copy
 
4
// of this software and associated documentation files (the "Software"), to
 
5
// deal in the Software without restriction, including without limitation the
 
6
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
 
7
// sell copies of the Software, and to permit persons to whom the Software is
 
8
// furnished to do so, subject to the following conditions:
 
9
// 
 
10
// The above copyright notice and this permission notice shall be included in
 
11
// all copies or substantial portions of the Software.
 
12
// 
 
13
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 
14
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 
15
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 
16
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 
17
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 
18
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 
19
// IN THE SOFTWARE.
 
20
 
 
21
#include "libmv/base/vector.h"
 
22
#include "libmv/logging/logging.h"
 
23
#include "libmv/multiview/fundamental.h"
 
24
#include "libmv/multiview/projection.h"
 
25
#include "libmv/numeric/levenberg_marquardt.h"
 
26
#include "libmv/numeric/numeric.h"
 
27
#include "libmv/simple_pipeline/reconstruction.h"
 
28
#include "libmv/simple_pipeline/tracks.h"
 
29
 
 
30
namespace libmv {
 
31
namespace {
 
32
 
 
33
void CoordinatesForMarkersInImage(const vector<Marker> &markers,
 
34
                                  int image,
 
35
                                  Mat *coordinates) {
 
36
  vector<Vec2> coords;
 
37
  for (int i = 0; i < markers.size(); ++i) {
 
38
    const Marker &marker = markers[i];
 
39
    if (markers[i].image == image) {
 
40
      coords.push_back(Vec2(marker.x, marker.y));
 
41
    }
 
42
  }
 
43
  coordinates->resize(2, coords.size());
 
44
  for (int i = 0; i < coords.size(); i++) {
 
45
    coordinates->col(i) = coords[i];
 
46
  }
 
47
}
 
48
 
 
49
void GetImagesInMarkers(const vector<Marker> &markers,
 
50
                        int *image1, int *image2) {
 
51
  if (markers.size() < 2) {
 
52
    return;
 
53
  }
 
54
  *image1 = markers[0].image;
 
55
  for (int i = 1; i < markers.size(); ++i) {
 
56
    if (markers[i].image != *image1) {
 
57
      *image2 = markers[i].image;
 
58
      return;
 
59
    }
 
60
  }
 
61
  *image2 = -1;
 
62
  LOG(FATAL) << "Only one image in the markers.";
 
63
}
 
64
 
 
65
}  // namespace
 
66
 
 
67
bool EuclideanReconstructTwoFrames(const vector<Marker> &markers,
 
68
                                   EuclideanReconstruction *reconstruction) {
 
69
  if (markers.size() < 16) {
 
70
    LG << "Not enough markers to initialize from two frames: " << markers.size();
 
71
    return false;
 
72
  }
 
73
 
 
74
  int image1, image2;
 
75
  GetImagesInMarkers(markers, &image1, &image2);
 
76
 
 
77
  Mat x1, x2;
 
78
  CoordinatesForMarkersInImage(markers, image1, &x1);
 
79
  CoordinatesForMarkersInImage(markers, image2, &x2);
 
80
 
 
81
  Mat3 F;
 
82
  NormalizedEightPointSolver(x1, x2, &F);
 
83
 
 
84
  // The F matrix should be an E matrix, but squash it just to be sure.
 
85
  Eigen::JacobiSVD<Mat3> svd(F, Eigen::ComputeFullU | Eigen::ComputeFullV);
 
86
 
 
87
  // See Hartley & Zisserman page 294, result 11.1, which shows how to get the
 
88
  // closest essential matrix to a matrix that is "almost" an essential matrix.
 
89
  double a = svd.singularValues()(0);
 
90
  double b = svd.singularValues()(1);
 
91
  double s = (a + b) / 2.0;
 
92
  LG << "Initial reconstruction's rotation is non-euclidean by "
 
93
     << (((a - b) / std::max(a, b)) * 100) << "%; singular values:"
 
94
     << svd.singularValues().transpose();
 
95
 
 
96
  Vec3 diag;
 
97
  diag << s, s, 0;
 
98
  Mat3 E = svd.matrixU() * diag.asDiagonal() * svd.matrixV().transpose();
 
99
 
 
100
  // Recover motion between the two images. Since this function assumes a
 
101
  // calibrated camera, use the identity for K.
 
102
  Mat3 R;
 
103
  Vec3 t;
 
104
  Mat3 K = Mat3::Identity();
 
105
  if (!MotionFromEssentialAndCorrespondence(E,
 
106
                                            K, x1.col(0),
 
107
                                            K, x2.col(0),
 
108
                                            &R, &t)) {
 
109
    LG << "Failed to compute R and t from E and K.";
 
110
    return false;
 
111
  }
 
112
 
 
113
  // Image 1 gets the reference frame, image 2 gets the relative motion.
 
114
  reconstruction->InsertCamera(image1, Mat3::Identity(), Vec3::Zero());
 
115
  reconstruction->InsertCamera(image2, R, t);
 
116
 
 
117
  LG << "From two frame reconstruction got:\nR:\n" << R
 
118
     << "\nt:" << t.transpose();
 
119
  return true;
 
120
}
 
121
 
 
122
namespace {
 
123
 
 
124
Mat3 DecodeF(const Vec9 &encoded_F) {
 
125
  // Decode F and force it to be rank 2.
 
126
  Map<const Mat3> full_rank_F(encoded_F.data(), 3, 3);
 
127
  Eigen::JacobiSVD<Mat3> svd(full_rank_F, Eigen::ComputeFullU | Eigen::ComputeFullV);
 
128
  Vec3 diagonal = svd.singularValues();
 
129
  diagonal(2) = 0;
 
130
  Mat3 F = svd.matrixU() * diagonal.asDiagonal() * svd.matrixV().transpose();
 
131
  return F;
 
132
}
 
133
 
 
134
// This is the stupidest way to refine F known to mankind, since it requires
 
135
// doing a full SVD of F at each iteration. This uses sampson error.
 
136
struct FundamentalSampsonCostFunction {
 
137
 public:
 
138
  typedef Vec   FMatrixType;
 
139
  typedef Vec9 XMatrixType;
 
140
 
 
141
  // Assumes markers are ordered by track.
 
142
  FundamentalSampsonCostFunction(const vector<Marker> &markers)
 
143
    : markers(markers) {}
 
144
 
 
145
  Vec operator()(const Vec9 &encoded_F) const {
 
146
    // Decode F and force it to be rank 2.
 
147
    Mat3 F = DecodeF(encoded_F);
 
148
    
 
149
    Vec residuals(markers.size() / 2);
 
150
    residuals.setZero();
 
151
    for (int i = 0; i < markers.size() / 2; ++i) {
 
152
      const Marker &marker1 = markers[2*i + 0];
 
153
      const Marker &marker2 = markers[2*i + 1];
 
154
      CHECK_EQ(marker1.track, marker2.track);
 
155
      Vec2 x1(marker1.x, marker1.y);
 
156
      Vec2 x2(marker2.x, marker2.y);
 
157
 
 
158
      residuals[i] = SampsonDistance(F, x1, x2);
 
159
    }
 
160
    return residuals;
 
161
  }
 
162
  const vector<Marker> &markers;
 
163
};
 
164
 
 
165
}  // namespace
 
166
 
 
167
bool ProjectiveReconstructTwoFrames(const vector<Marker> &markers,
 
168
                                    ProjectiveReconstruction *reconstruction) {
 
169
  if (markers.size() < 16) {
 
170
    return false;
 
171
  }
 
172
 
 
173
  int image1, image2;
 
174
  GetImagesInMarkers(markers, &image1, &image2);
 
175
 
 
176
  Mat x1, x2;
 
177
  CoordinatesForMarkersInImage(markers, image1, &x1);
 
178
  CoordinatesForMarkersInImage(markers, image2, &x2);
 
179
 
 
180
  Mat3 F;
 
181
  NormalizedEightPointSolver(x1, x2, &F);
 
182
 
 
183
  // XXX Verify sampson distance.
 
184
#if 0
 
185
  // Refine the resulting projection fundamental matrix using Sampson's
 
186
  // approximation of geometric error. This avoids having to do a full bundle
 
187
  // at the cost of some accuracy.
 
188
  //
 
189
  // TODO(keir): After switching to a better bundling library, use a proper
 
190
  // full bundle adjust here instead of this lame bundle adjustment.
 
191
  typedef LevenbergMarquardt<FundamentalSampsonCostFunction> Solver;
 
192
 
 
193
  FundamentalSampsonCostFunction fundamental_cost(markers);
 
194
 
 
195
  // Pack the initial P matrix into a size-12 vector..
 
196
  Vec9 encoded_F = Map<Vec9>(F.data(), 3, 3);
 
197
 
 
198
  Solver solver(fundamental_cost);
 
199
 
 
200
  Solver::SolverParameters params;
 
201
  Solver::Results results = solver.minimize(params, &encoded_F);
 
202
  // TODO(keir): Check results to ensure clean termination.
 
203
 
 
204
  // Recover F from the minimization.
 
205
  F = DecodeF(encoded_F);
 
206
#endif
 
207
 
 
208
  // Image 1 gets P = [I|0], image 2 gets arbitrary P.
 
209
  Mat34 P1 = Mat34::Zero();
 
210
  P1.block<3, 3>(0, 0) = Mat3::Identity();
 
211
  Mat34 P2;
 
212
  ProjectionsFromFundamental(F, &P1, &P2);
 
213
 
 
214
  reconstruction->InsertCamera(image1, P1);
 
215
  reconstruction->InsertCamera(image2, P2);
 
216
 
 
217
  LG << "From two frame reconstruction got P2:\n" << P2;
 
218
  return true;
 
219
}
 
220
}  // namespace libmv