griefith/test
1
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases Wiki Activity

Update libmv to current upstream version

Browse Source 
- Solves some strict compilation warning
- Style/code cleanup
This commit is contained in:
Sergey Sharybin
2013-04-05 09:23:20 +00:00
parent 43b61fb8bd
commit 898ba93a12
77 changed files with 1066 additions and 1039 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
extern/libmv/CMakeLists.txt vendored
View File

@@ -45,12 +45,12 @@ set(SRC
libmv/multiview/euclidean_resection.cc
libmv/multiview/fundamental.cc
libmv/multiview/homography.cc
libmv/multiview/panography.cc
libmv/multiview/projection.cc
libmv/multiview/triangulation.cc
libmv/numeric/numeric.cc
libmv/numeric/poly.cc
libmv/simple_pipeline/bundle.cc
libmv/simple_pipeline/callbacks.cc
libmv/simple_pipeline/camera_intrinsics.cc
libmv/simple_pipeline/detect.cc
libmv/simple_pipeline/initialize_reconstruction.cc

431
extern/libmv/ChangeLog vendored
View File

@@ -1,3 +1,169 @@
commit e3b2bccba6145290738a6677c14f7369ec7a38cd
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Thu Apr 4 02:59:58 2013 +0600
Suppress strict compiler warnings in glags/glog libraries
commit 5fca459adcf0a3419fa9cd8d983dc2c952d02647
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Thu Apr 4 01:20:18 2013 +0600
Lint cleanup, mostly white space and line width.
Also moved own includes to the top of files.
Should be no functional changes :)
commit 9a9dd458a622928b91dbd3c79900577923283838
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Fri Mar 29 00:20:29 2013 +0600
Fix for TransformTracks in uncalibrated pipeline
Transformation matrix was completely ignored by
TransformTracks() and final marker coordinate
exactly matched it's source coordinates.
Seems to be just a typo in vector usage: need to
use "b" (which is transformed one) instead of "a"
when converting projective coordinate to 2D space.
commit d35766cc9901609e32f4d80faba715695bea3c40
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Fri Mar 29 00:19:11 2013 +0600
Revert part of e2eb58c4230f94ef0c72fb4005e5434088d52e80
That commit included one change which shall have been
go as separate commit with more detailed description.
commit e8d71b4e96fd78eb60773b6557d66da672e65753
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Wed Mar 27 20:37:05 2013 +0600
Silenced more warnings
- Added includes of own header to fast implementation files.
- Camera intrinsics wouldn't complain about unknown pragma when
building without OpenMP support.
TODO: Make it a CMake option to build libmv with OpenMP support.
Currently multi-threaded intrinsics only available when
using custom CMake rules for bundled libmv version
(as it's done in Blender).
commit ad442812654f270dc088394410fda1b81b8dc450
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Wed Mar 27 20:18:51 2013 +0600
Multithreaded camera intrinsics
Implemented multithreaded buffer (un)distortion
for camera intrinsics using OpenMP.
By default, (un)distortion is single-threaded,
but it is possible to as CameraIntrinsics to
use more threads by calling SetThreads method.
commit c88b4881096174a16a9f9e6fc2c9dcad3e255b25
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Wed Mar 27 18:45:09 2013 +0600
Movie functions implementation from panography header
into own CC implementation file.
Before this all panography functions were declared as
static, which is not so much useful from re-useability
point of view.
commit 2d2faf9104bc035722cff6775e1b8e7c93143aba
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Wed Mar 27 18:37:36 2013 +0600
Build shared Ceres library only if BUILD_SHARED_LIBS is enabled
commit daa3ddd3260ccaf2bf9c72eadb89213d91e549ec
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Wed Mar 27 18:21:52 2013 +0600
Update Ceres to upstream version 1.5.0
commit cf5dc678878345ea3f221ce50cb2b9e539c2ab38
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Wed Mar 27 15:06:24 2013 +0600
Code cleanup: removed more deprecated FFmpeg API usage
This time in qt-tracker application.
commit e2eb58c4230f94ef0c72fb4005e5434088d52e80
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Tue Mar 26 17:19:51 2013 +0600
Code cleanup: silent unused variables warnings
commit af89bb24667e39b7e655173ea807fdcfbeef4422
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Tue Mar 26 16:54:14 2013 +0600
Code cleanup: no need to declare empty body for ProgressUpdateCallback:invoke
Make force this method to be overridden by derivative classes.
Also removed currently unneeded callbacks.cc.
commit 0441d4ee06fad0219256a5704f931eec916a3868
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Tue Mar 26 16:37:27 2013 +0600
Code cleanup: silent type narrowing in qt-tracker
commit cd4b61c976448d0fdedefb3ed4b21d70e078f94b
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Tue Mar 26 16:26:39 2013 +0600
Changes to unit testing
- Move ceres test binaries to ${LIBMV_TESTS_OUTPUT_DIR}/ceres,
so they don't mess with libmv's application binaries and
tests.
- Removed ceres_ prefix from ceres unit tests, only use this
prefix for targets (targets need to be unique name).
- Added unit tests data for ceres, otherwise system_test fails.
- Restored "test" makefile target.
commit cf704ada08acc8b26167e7bfb3e1e88fd278de23
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Tue Mar 26 15:01:15 2013 +0600
Code cleanup: use rw-rw-r-- mode for source files
commit 64b31e3e43acb52aaf6f591b9d1c2449bf6ef3bd
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Tue Mar 26 14:57:46 2013 +0600
Code cleanup: don't use deprecated FFmpeg API functions
commit 2a3676499548ad5dba5a5c5eadf3bb71e640b612
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Tue Mar 5 17:40:52 2013 +0600
Switch from DENSE_NORMAL_CHOLESKY to DENSE_QR
DENSE_QR is better behaved numerically and after recent
changes from Sameer there's no big difference in speed.
commit bcb920df02133da5b7e55fbc74edb9222004eecc
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Tue Mar 5 17:15:43 2013 +0600
Update Ceres to 1.5RC3
It brings optimization of DENSE_QR and DENSE_SCHUR solvers.
commit 473996468a4e67e7c860169181a4ff31ce9b8c80
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Fri Mar 1 17:44:54 2013 +0600
@@ -533,268 +699,3 @@ Date: Fri Jun 8 17:42:17 2012 +0000
aborts immediately. The workaround for now is to disable the
correlation checking, and examine your tracks carefully. A
fix will get added shortly.
commit 81d028f13738ebe2304287dfce90e91bc782e2cf
Author: Keir Mierle <mierle@gmail.com>
Date: Fri May 18 20:04:43 2012 +0000
Remove an unnecessary template<> line in libmv. Convert debug logs to LG.
commit 238aaba241ef99995d254aadc974db719da04b96
Author: Keir Mierle <mierle@gmail.com>
Date: Fri May 18 12:05:10 2012 +0000
Support normalization in the tracking prepass
The last tracker commit added normalized tracking. This makes
tracking patches undergoing uniform illumination change easier.
However, the prepass which computes a quick translation-only
estimate of the warp did not take this into account. This commit
fixes that.
This works reasonably well but in some examples the brute
initialization fails. I suspect this is due to the warped template
estimate in the current frame being too different from the
original, so there are multiple peaks in the normalized-SAD
correlation function.
The solution is to use the previous frame for the brute
initialization and the keyframe for refinement, but that requires
architecture changes.
commit 981ca4f6a679cd9ac3d086eae3cd946ce72ca8a5
Author: Keir Mierle <mierle@gmail.com>
Date: Fri May 18 02:12:47 2012 +0000
Add light-normalized tracking to the planar tracker
This commit adds the ability to normalize patterns by their
average value while tracking, to make them invariant to global
illumination changes.
To see this in action, check out the "Lobby" scene from Hollywood
VFX. If you track the markers that are shadowed by the actress,
previously they would not track. With the scale adaption on, the
tracker would shrink the area to compensate for the changed
illumination, losing the track. With "Normalize" turned on, the
patch is correctly tracked and scale is maintained.
A remaining problem is that only the Ceres cost function is
updated to handle the normalization. The brute translation search
does not take this into account. Perhaps "Prepass" (see below)
should get disabled if normalization is enabled until I fix the
prepass to normalize as well.
There are a few other changes:
- Bail out of the sampling loop early if the mask is zero; this
saves expensive samples of the image derivatives.
- Fix a bug where the mask was ignored when sampling in the cost
functor.
commit e9384b15fb2a6a5b81346d5758fa136f0911e945
Author: Keir Mierle <mierle@gmail.com>
Date: Thu May 17 23:53:32 2012 +0000
Implement support for affine tracking in the planar tracker; cleanups.
commit 021d41eed8b4ce6a4e37786ccd357ed5dc83a13f
Author: Keir Mierle <mierle@gmail.com>
Date: Thu May 17 21:26:06 2012 +0000
For the planar tracker, initialize the warp from the four correspondences
after brute force translation search.
commit 003d1bf6145cfd30938b35f6e10d43708dbf916c
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Thu Dec 6 16:56:01 2012 +0600
Correction to region tracker options initialization.
Based on patch from Keir to Blender:
https://svn.blender.org/svnroot/bf-blender/branches/soc-2011-tomato@46743
commit 6af47b218cfdf5219f0ebb3cb95459817cf9abf2
Author: Keir Mierle <mierle@gmail.com>
Date: Thu May 17 02:31:52 2012 +0000
Add new planar tracker features and use the new planar API
This commit removes the use of the legacy RegionTracker API from
Blender, and replaces it with the new TrackRegion API. This also
adds several features to the planar tracker in libmv:
- Do a brute-force initialization of tracking similar to "Hybrid"
mode in the stable release, but using all floats. This is slower
but more accurate. It is still necessary to evaluate if the
performance loss is worth it. In particular, this change is
necessary to support high bit depth imagery.
- Add support for masks over the search window. This is a step
towards supporting user-defined tracker masks. The tracker masks
will make it easy for users to make a mask for e.g. a ball.
- Add Pearson product moment correlation coefficient checking (aka
"Correlation" in the UI. This causes tracking failure if the
tracked patch is not linearly related to the template.
- Add support for warping a few points in addition to the supplied
points. This is useful because the tracking code deliberately
does not expose the underlying warp representation. Instead,
warps are specified in an aparametric way via the correspondences.
- Remove the "num_samples_xy" concept and replace it with
automatic determination of the number of samples. This makes the
API easier for users.
- Fix various bugs in the parameterizations.
There remains a bug with subpixel precision tracking when in
"keyframe" mode; this will get fixed shortly.
commit 16a46db104468cec80bd31ca9d5f8bffbe3e003e
Author: Keir Mierle <mierle@gmail.com>
Date: Mon May 14 12:15:38 2012 +0000
"Efficient Second-order Minimization" for the planar tracker
This implements the "Efficient Second-order Minimization"
scheme, as supported by the existing translation tracker.
This increases the amount of per-iteration work, but
decreases the number of iterations required to converge and
also increases the size of the basin of attraction for the
optimization.
commit 23243b1b1f3e1ab3ef862b47bca06ee876ac2cf4
Author: Keir Mierle <mierle@gmail.com>
Date: Sun May 13 23:08:56 2012 +0000
Add a planar tracking implementation to libmv
This adds a new planar tracking implementation to libmv. The
tracker is based on Ceres[1], the new nonlinear minimizer that
myself and Sameer released from Google as open source. Since
the motion model is more involved, the interface is
different than the RegionTracker interface used previously
in Blender.
The ESM tracker, also known as the KLT tracker in the UI, is
temporarily changed to use the new Ceres-based planar
tracker in translation-only mode. Currently it is a bit
slower than ESM and also doesn't have all the bells and
whistles implemented. Those will come soon. Longer term,
both trackers will remain since Ceres is unlikely to be as
fast as ESM for pure translation solving, due to its
generality.
The next step is to implement a new tracking UI. The current
UI assumes a translational motion model; the new one must
support arbitrary perspective transforms of the pattern
regions.
[1] http://code.google.com/p/ceres-solver
commit 52be92b53eb4decb1a316690b162196f227cc441
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Thu Dec 6 16:06:08 2012 +0600
Initial Ceres integration
Currently only put sources to src/third_party/ceres and made sure they're
not giving compilation issues.
Used Ceres upstream version 1.3.0.
Needed to make some modifications to it's CMakeLists.txt also to glog and
fglags. They're described in README.libmv of this libraries.
Basically:
- Added -fPIC to glog/gflags, so shared ceres library could be linked
statically against this libraries.
- Tweaked Ceres's build rules to use needed libraries from libmv's
third_party folder.
commit b13f9d13122e091cb85855c2094386ccdef6e5a4
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Wed Dec 5 19:05:34 2012 +0600
Update Eigen to version 3.1.2
Mainly because of lots of warnings generating by gcc-4.7 which are
resolved in newer eigen version.
commit 1f0dd94e8e37d3fe2df89282ec16a6a685fdde0b
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Fri May 25 16:36:44 2012 +0600
- Added avutil to qt-tracker linking when building with FFmpeg support.
On some platforms it seems to be required
- Synchronized QT Creator project for qt-tracker with changes in sources,
so no it might be compiled from QT Creator.
commit b813dbe3f46bbbc7e73ac791d4665622e4fc7ba5
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Wed May 9 19:01:10 2012 +0600
Modal solver: Detect rigid transformation between initial frame and current
instead of detecting it between two neighbour frames.
This prevents accumulation of error and seems to be working better in footages i've tested.
commit 9254621c76daaf239ec1f535e197ca792eea97b6
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Wed May 9 18:57:00 2012 +0600
Backport changes made by Keir in Blender:
- Enhance logging in libmv's trackers.
- Cleanups in brute_region_tracker.cc.
commit d9c56b9d3c63f886d83129ca0ebed1e76d9c93d7
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Fri Apr 27 16:20:41 2012 +0600
Fixes for MinGW64 support by Caleb Joseph with slight modifications by Antony Riakiotakis
- Functions snprintf and sincos shouldn't be redefined for MinGW64
- Type pid_t shouldn't be re-defined for MinGW64
commit e1902b6938676011607ac99986b8b140bdbf090e
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Fri Apr 27 16:04:19 2012 +0600
Fixes for Qt calibration tool
- Passing directory with images via command line argument now isn't
required -- it there's no such directory specified standard open
dialog might be used for this (before application used to abort
due to accessing to non-existing list element).
- Conversion of source images to grayscale now happens correct.
It was needed to build grayscale palette for 8bit indexed buffer.
commit 05f1a0a78ad8ff6646d1e8da97e6f7575b891536
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Sat Apr 14 17:21:29 2012 +0600
Make QtTracker compilable again porting it to recent API change and code cleanup:
- It was using SAD tracker with own API, now it's using standard RegionTracker API
which should make it easier to switch between different trackers.
- Restored LaplaceFilter from old SAD module which convolves images with the
discrete laplacian operator.
commit a44312a7beb2963b8e3bf8015c516d2eff40cc3d
Author: Sergey Sharybin <sergey.vfx@gmail.com>
Date: Thu Apr 12 13:56:02 2012 +0600
Added solver for modal camera motion, currently supports only tripod solving
This solver is intended to deal with such camera motions as tripod and panning,
where it's impossible to reconstruct exact position of markers in 3d view.
It projects markers onto sphere and uses rigid registration of rotation to
find rotation angles which makes bundles from previous and current frame be
as closest as it's possible.

2
extern/libmv/files.txt vendored
View File

@@ -21,6 +21,7 @@ libmv/multiview/homography.cc
libmv/multiview/homography.h
libmv/multiview/homography_parameterization.h
libmv/multiview/nviewtriangulation.h
libmv/multiview/panography.cc
libmv/multiview/panography.h
libmv/multiview/projection.cc
libmv/multiview/projection.h
@@ -36,7 +37,6 @@ libmv/numeric/poly.cc
libmv/numeric/poly.h
libmv/simple_pipeline/bundle.cc
libmv/simple_pipeline/bundle.h
libmv/simple_pipeline/callbacks.cc
libmv/simple_pipeline/callbacks.h
libmv/simple_pipeline/camera_intrinsics.cc
libmv/simple_pipeline/camera_intrinsics.h

2
extern/libmv/libmv/base/vector.h vendored
View File

@@ -39,7 +39,7 @@ namespace libmv {
// - doesn't support iterators.
// - impede compatibility with code using STL.
// - the STL already provide support for custom allocators
// it could be replaced with a simple
// it could be replaced with a simple
// template <T> class vector : std::vector<T, aligned_allocator> {} declaration
// provided it doesn't break code relying on libmv::vector specific behavior
template <typename T,

2
extern/libmv/libmv/base/vector_utils.h vendored
View File

@@ -31,4 +31,4 @@ void DeleteElements(Array *array) {
array->clear();
}
#endif // LIBMV_BASE_VECTOR_UTILS_H_
#endif // LIBMV_BASE_VECTOR_UTILS_H_

10
extern/libmv/libmv/image/array_nd.cc vendored
View File

@@ -35,8 +35,8 @@ void FloatArrayToScaledByteArray(const Array3Df &float_array,
for (int i = 0; i < float_array.Height(); ++i) {
for (int j = 0; j < float_array.Width(); ++j) {
for (int k = 0; k < float_array.Depth(); ++k) {
minval = std::min(minval, float_array(i,j,k));
maxval = std::max(maxval, float_array(i,j,k));
minval = std::min(minval, float_array(i, j, k));
maxval = std::max(maxval, float_array(i, j, k));
}
}
}
@@ -47,8 +47,8 @@ void FloatArrayToScaledByteArray(const Array3Df &float_array,
for (int i = 0; i < float_array.Height(); ++i) {
for (int j = 0; j < float_array.Width(); ++j) {
for (int k = 0; k < float_array.Depth(); ++k) {
float unscaled = (float_array(i,j,k) - minval) / (maxval - minval);
(*byte_array)(i,j,k) = (unsigned char)(255 * unscaled);
float unscaled = (float_array(i, j, k) - minval) / (maxval - minval);
(*byte_array)(i, j, k) = (unsigned char)(255 * unscaled);
}
}
}
@@ -60,7 +60,7 @@ void ByteArrayToScaledFloatArray(const Array3Du &byte_array,
for (int i = 0; i < byte_array.Height(); ++i) {
for (int j = 0; j < byte_array.Width(); ++j) {
for (int k = 0; k < byte_array.Depth(); ++k) {
(*float_array)(i,j,k) = float(byte_array(i,j,k)) / 255.0f;
(*float_array)(i, j, k) = float(byte_array(i, j, k)) / 255.0f;
}
}
}

44
extern/libmv/libmv/image/array_nd.h vendored
View File

@@ -57,9 +57,13 @@ class ArrayND : public BaseArray {
ArrayND(int s0) : data_(NULL), own_data(true) { Resize(s0); }
ArrayND(int s0, int s1) : data_(NULL), own_data(true) { Resize(s0, s1); }
ArrayND(int s0, int s1, int s2) : data_(NULL), own_data(true) { Resize(s0, s1, s2); }
ArrayND(int s0, int s1, int s2) : data_(NULL), own_data(true) {
Resize(s0, s1, s2);
}
ArrayND(T* data, int s0, int s1, int s2) : data_(data), own_data(false) { Resize(s0, s1, s2); }
ArrayND(T* data, int s0, int s1, int s2) : data_(data), own_data(false) {
Resize(s0, s1, s2);
}
/// Destructor deletes pixel data.
~ArrayND() {
@@ -93,7 +97,7 @@ class ArrayND : public BaseArray {
for (int i = N - 1; i > 0; --i) {
strides_(i - 1) = strides_(i) * shape_(i);
}
if(own_data) {
if (own_data) {
delete [] data_;
data_ = NULL;
if (Size() > 0) {
@@ -103,7 +107,7 @@ class ArrayND : public BaseArray {
}
template<typename D>
void ResizeLike(const ArrayND<D,N> &other) {
void ResizeLike(const ArrayND<D, N> &other) {
Resize(other.Shape());
}
@@ -136,12 +140,12 @@ class ArrayND : public BaseArray {
/// Resize a 3D array to shape (s0,s1,s2).
void Resize(int s0, int s1, int s2) {
assert(N == 3);
int shape[] = {s0,s1,s2};
int shape[] = {s0, s1, s2};
Resize(shape);
}
template<typename D>
void CopyFrom(const ArrayND<D,N> &other) {
void CopyFrom(const ArrayND<D, N> &other) {
ResizeLike(other);
T *data = Data();
const D *other_data = other.Data();
@@ -238,7 +242,7 @@ class ArrayND : public BaseArray {
T &operator()(int i0, int i1) {
assert(0 <= i0 && i0 < Shape(0));
assert(0 <= i1 && i1 < Shape(1));
return *( Data() + Offset(i0,i1) );
return *(Data() + Offset(i0, i1));
}
/// 3D specialization.
@@ -246,29 +250,29 @@ class ArrayND : public BaseArray {
assert(0 <= i0 && i0 < Shape(0));
assert(0 <= i1 && i1 < Shape(1));
assert(0 <= i2 && i2 < Shape(2));
return *( Data() + Offset(i0,i1,i2) );
return *(Data() + Offset(i0, i1, i2));
}
/// Return a constant reference to the element at position index.
const T &operator()(const Index &index) const {
return *( Data() + Offset(index) );
return *(Data() + Offset(index));
}
/// 1D specialization.
const T &operator()(int i0) const {
return *( Data() + Offset(i0) );
return *(Data() + Offset(i0));
}
/// 2D specialization.
const T &operator()(int i0, int i1) const {
assert(0 <= i0 && i0 < Shape(0));
assert(0 <= i1 && i1 < Shape(1));
return *( Data() + Offset(i0,i1) );
return *(Data() + Offset(i0, i1));
}
/// 3D specialization.
const T &operator()(int i0, int i1, int i2) const {
return *( Data() + Offset(i0,i1,i2) );
return *(Data() + Offset(i0, i1, i2));
}
/// True if index is inside array.
@@ -343,14 +347,14 @@ class Array3D : public ArrayND<T, 3> {
Array3D()
: Base() {
}
Array3D(int height, int width, int depth=1)
Array3D(int height, int width, int depth = 1)
: Base(height, width, depth) {
}
Array3D(T* data, int height, int width, int depth=1)
Array3D(T* data, int height, int width, int depth = 1)
: Base(data, height, width, depth) {
}
void Resize(int height, int width, int depth=1) {
void Resize(int height, int width, int depth = 1) {
Base::Resize(height, width, depth);
}
@@ -376,12 +380,12 @@ class Array3D : public ArrayND<T, 3> {
T &operator()(int i0, int i1, int i2 = 0) {
assert(0 <= i0 && i0 < Height());
assert(0 <= i1 && i1 < Width());
return Base::operator()(i0,i1,i2);
return Base::operator()(i0, i1, i2);
}
const T &operator()(int i0, int i1, int i2 = 0) const {
assert(0 <= i0 && i0 < Height());
assert(0 <= i1 && i1 < Width());
return Base::operator()(i0,i1,i2);
return Base::operator()(i0, i1, i2);
}
};
@@ -413,13 +417,13 @@ void ByteArrayToScaledFloatArray(const Array3Du &byte_array,
Array3Df *float_array);
template <typename AArrayType, typename BArrayType, typename CArrayType>
void MultiplyElements( const AArrayType &a,
void MultiplyElements(const AArrayType &a,
const BArrayType &b,
CArrayType *c ) {
CArrayType *c) {
// This function does an element-wise multiply between
// the two Arrays A and B, and stores the result in C.
// A and B must have the same dimensions.
assert( a.Shape() == b.Shape() );
assert(a.Shape() == b.Shape());
c->ResizeLike(a);
// To perform the multiplcation, a "current" index into the N-dimensions of

60
extern/libmv/libmv/image/convolve.cc vendored
View File

@@ -18,10 +18,11 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/image/convolve.h"
#include <cmath>
#include "libmv/image/image.h"
#include "libmv/image/convolve.h"
namespace libmv {
@@ -118,7 +119,7 @@ void Convolve(const Array3Df &in,
// fast path.
int half_width = kernel.size() / 2;
switch (half_width) {
#define static_convolution( size ) case size: \
#define static_convolution(size) case size: \
FastConvolve<size, vertical>(kernel, width, height, src, src_stride, \
src_line_stride, dst, dst_stride); break;
static_convolution(1)
@@ -136,13 +137,15 @@ void Convolve(const Array3Df &in,
double sum = 0;
// Slow path: this loop cannot be unrolled.
for (int k = -dynamic_size; k <= dynamic_size; ++k) {
if(vertical) {
if (vertical) {
if (y + k >= 0 && y + k < height) {
sum += src[k * src_line_stride] * kernel(2 * dynamic_size - (k + dynamic_size));
sum += src[k * src_line_stride] *
kernel(2 * dynamic_size - (k + dynamic_size));
}
} else {
if (x + k >= 0 && x + k < width) {
sum += src[k * src_stride] * kernel(2 * dynamic_size - (k + dynamic_size));
sum += src[k * src_stride] *
kernel(2 * dynamic_size - (k + dynamic_size));
}
}
}
@@ -234,25 +237,25 @@ void BoxFilterHorizontal(const Array3Df &in,
int half_width = (window_size - 1) / 2;
for (int k = 0; k < in.Depth(); ++k) {
for (int i=0; i<in.Height(); ++i) {
for (int i = 0; i < in.Height(); ++i) {
float sum = 0;
// Init sum.
for (int j=0; j<half_width; ++j) {
for (int j = 0; j < half_width; ++j) {
sum += in(i, j, k);
}
// Fill left border.
for (int j=0; j < half_width + 1; ++j) {
for (int j = 0; j < half_width + 1; ++j) {
sum += in(i, j + half_width, k);
out(i, j, k) = sum;
}
// Fill interior.
for (int j = half_width + 1; j<in.Width()-half_width; ++j) {
for (int j = half_width + 1; j < in.Width()-half_width; ++j) {
sum -= in(i, j - half_width - 1, k);
sum += in(i, j + half_width, k);
out(i, j, k) = sum;
}
// Fill right border.
for (int j = in.Width() - half_width; j<in.Width(); ++j) {
for (int j = in.Width() - half_width; j < in.Width(); ++j) {
sum -= in(i, j - half_width - 1, k);
out(i, j, k) = sum;
}
@@ -271,22 +274,22 @@ void BoxFilterVertical(const Array3Df &in,
for (int j = 0; j < in.Width(); ++j) {
float sum = 0;
// Init sum.
for (int i=0; i<half_width; ++i) {
for (int i = 0; i < half_width; ++i) {
sum += in(i, j, k);
}
// Fill left border.
for (int i=0; i < half_width + 1; ++i) {
for (int i = 0; i < half_width + 1; ++i) {
sum += in(i + half_width, j, k);
out(i, j, k) = sum;
}
// Fill interior.
for (int i = half_width + 1; i<in.Height()-half_width; ++i) {
for (int i = half_width + 1; i < in.Height()-half_width; ++i) {
sum -= in(i - half_width - 1, j, k);
sum += in(i + half_width, j, k);
out(i, j, k) = sum;
}
// Fill right border.
for (int i = in.Height() - half_width; i<in.Height(); ++i) {
for (int i = in.Height() - half_width; i < in.Height(); ++i) {
sum -= in(i - half_width - 1, j, k);
out(i, j, k) = sum;
}
@@ -302,18 +305,23 @@ void BoxFilter(const Array3Df &in,
BoxFilterVertical(tmp, box_width, out);
}
void LaplaceFilter(unsigned char* src, unsigned char* dst, int width, int height, int strength) {
for(int y=1; y<height-1; y++) for(int x=1; x<width-1; x++) {
const unsigned char* s = &src[y*width+x];
int l = 128 +
s[-width-1] + s[-width] + s[-width+1] +
s[1] - 8*s[0] + s[1] +
s[ width-1] + s[ width] + s[ width+1] ;
int d = ((256-strength)*s[0] + strength*l) / 256;
if(d < 0) d=0;
if(d > 255) d=255;
dst[y*width+x] = d;
}
void LaplaceFilter(unsigned char* src,
unsigned char* dst,
int width,
int height,
int strength) {
for (int y = 1; y < height-1; y++)
for (int x = 1; x < width-1; x++) {
const unsigned char* s = &src[y*width+x];
int l = 128 +
s[-width-1] + s[-width] + s[-width+1] +
s[1] - 8*s[0] + s[1] +
s[ width-1] + s[ width] + s[ width+1];
int d = ((256-strength)*s[0] + strength*l) / 256;
if (d < 0) d=0;
if (d > 255) d=255;
dst[y*width+x] = d;
}
}
} // namespace libmv

6
extern/libmv/libmv/image/convolve.h vendored
View File

@@ -95,7 +95,11 @@ void BoxFilter(const FloatImage &in,
\note Make sure the search region is filtered with the same strength as the pattern.
*/
void LaplaceFilter(unsigned char* src, unsigned char* dst, int width, int height, int strength);
void LaplaceFilter(unsigned char* src,
unsigned char* dst,
int width,
int height,
int strength);
} // namespace libmv

7
extern/libmv/libmv/image/correlation.h vendored
View File

@@ -26,9 +26,10 @@
namespace libmv {
inline double PearsonProductMomentCorrelation(Array3Df image_and_gradient1_sampled,
Array3Df image_and_gradient2_sampled,
int width) {
inline double PearsonProductMomentCorrelation(
Array3Df image_and_gradient1_sampled,
Array3Df image_and_gradient2_sampled,
int width) {
double sX = 0, sY = 0, sXX = 0, sYY = 0, sXY = 0;
for (int r = 0; r < width; ++r) {

17
extern/libmv/libmv/image/image.h vendored
View File

@@ -60,8 +60,7 @@ class Image {
// Size in bytes that the image takes in memory.
int MemorySizeInBytes() {
int size;
switch (array_type_)
{
switch (array_type_) {
case BYTE:
size = reinterpret_cast<Array3Du *>(array_)->MemorySizeInBytes();
break;
@@ -83,8 +82,7 @@ class Image {
}
~Image() {
switch (array_type_)
{
switch (array_type_) {
case BYTE:
delete reinterpret_cast<Array3Du *>(array_);
@@ -109,23 +107,22 @@ class Image {
Image& operator= (const Image& f) {
if (this != &f) {
array_type_ = f.array_type_;
switch (array_type_)
{
switch (array_type_) {
case BYTE:
delete reinterpret_cast<Array3Du *>(array_);
array_ = new Array3Du( *(Array3Du *)f.array_);
array_ = new Array3Du(*(Array3Du *)f.array_);
break;
case FLOAT:
delete reinterpret_cast<Array3Df *>(array_);
array_ = new Array3Df( *(Array3Df *)f.array_);
array_ = new Array3Df(*(Array3Df *)f.array_);
break;
case INT:
delete reinterpret_cast<Array3Di *>(array_);
array_ = new Array3Di( *(Array3Di *)f.array_);
array_ = new Array3Di(*(Array3Di *)f.array_);
break;
case SHORT:
delete reinterpret_cast<Array3Ds *>(array_);
array_ = new Array3Ds( *(Array3Ds *)f.array_);
array_ = new Array3Ds(*(Array3Ds *)f.array_);
break;
default:
assert(0);

8
extern/libmv/libmv/image/sample.h vendored
View File

@@ -71,8 +71,8 @@ inline T SampleLinear(const Array3D<T> &image, float y, float x, int v = 0) {
const T im21 = image(y2, x1, v);
const T im22 = image(y2, x2, v);
return T( dy * ( dx * im11 + (1.0 - dx) * im12 ) +
(1 - dy) * ( dx * im21 + (1.0 - dx) * im22 ));
return T( dy * (dx * im11 + (1.0 - dx) * im12) +
(1 - dy) * (dx * im21 + (1.0 - dx) * im22));
}
/// Linear interpolation, of all channels. The sample is assumed to have the
@@ -95,8 +95,8 @@ inline void SampleLinear(const Array3D<T> &image, float y, float x, T *sample) {
const T im21 = image(y2, x1, i);
const T im22 = image(y2, x2, i);
sample[i] = T( dy * ( dx * im11 + (1.0 - dx) * im12 ) +
(1 - dy) * ( dx * im21 + (1.0 - dx) * im22 ));
sample[i] = T( dy * (dx * im11 + (1.0 - dx) * im12) +
(1 - dy) * (dx * im21 + (1.0 - dx) * im22));
}
}

8
extern/libmv/libmv/image/tuple.h vendored
View File

@@ -37,10 +37,10 @@ class Tuple {
Tuple(D *values) { Reset(values); }
template <typename D>
Tuple(const Tuple<D,N> &b) { Reset(b); }
Tuple(const Tuple<D, N> &b) { Reset(b); }
template <typename D>
Tuple& operator=(const Tuple<D,N>& b) {
Tuple& operator=(const Tuple<D, N>& b) {
Reset(b);
return *this;
}
@@ -50,14 +50,14 @@ class Tuple {
template <typename D>
void Reset(D *values) {
for(int i=0;i<N;i++) {
for (int i = 0;i < N; i++) {
data_[i] = T(values[i]);
}
}
// Set all tuple values to the same thing.
void Reset(T value) {
for(int i=0;i<N;i++) {
for (int i = 0;i < N; i++) {
data_[i] = value;
}
}

6
extern/libmv/libmv/multiview/conditioning.cc vendored
View File

@@ -22,7 +22,7 @@
#include "libmv/multiview/projection.h"
namespace libmv {
// HZ 4.4.4 pag.109: Point conditioning (non isotropic)
void PreconditionerFromPoints(const Mat &points, Mat3 *T) {
Vec mean, variance;
@@ -66,7 +66,7 @@ void ApplyTransformationToPoints(const Mat &points,
const Mat3 &T,
Mat *transformed_points) {
int n = points.cols();
transformed_points->resize(2,n);
transformed_points->resize(2, n);
Mat3X p(3, n);
EuclideanToHomogeneous(points, &p);
p = T * p;
@@ -96,4 +96,4 @@ void UnnormalizerI::Unnormalize(const Mat3 &T1, const Mat3 &T2, Mat3 *H) {
*H = T2.inverse() * (*H) * T1;
}
} // namespace libmv
} // namespace libmv

5
extern/libmv/libmv/multiview/conditioning.h vendored
View File

@@ -54,7 +54,6 @@ struct UnnormalizerT {
static void Unnormalize(const Mat3 &T1, const Mat3 &T2, Mat3 *H);
};
} //namespace libmv
} // namespace libmv
#endif // LIBMV_MULTIVIEW_CONDITIONNING_H_
#endif // LIBMV_MULTIVIEW_CONDITIONNING_H_

126
extern/libmv/libmv/multiview/euclidean_resection.cc vendored
View File

@@ -18,6 +18,8 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/multiview/euclidean_resection.h"
#include <cmath>
#include <limits>
@@ -26,7 +28,6 @@
#include "libmv/base/vector.h"
#include "libmv/logging/logging.h"
#include "libmv/multiview/euclidean_resection.h"
#include "libmv/multiview/projection.h"
namespace libmv {
@@ -34,7 +35,7 @@ namespace euclidean_resection {
typedef unsigned int uint;
bool EuclideanResection(const Mat2X &x_camera,
bool EuclideanResection(const Mat2X &x_camera,
const Mat3X &X_world,
Mat3 *R, Vec3 *t,
ResectionMethod method,
@@ -52,7 +53,7 @@ bool EuclideanResection(const Mat2X &x_camera,
return false;
}
bool EuclideanResection(const Mat &x_image,
bool EuclideanResection(const Mat &x_image,
const Mat3X &X_world,
const Mat3 &K,
Mat3 *R, Vec3 *t,
@@ -80,11 +81,11 @@ void AbsoluteOrientation(const Mat3X &X,
// Normalize the two point sets.
Mat3X Xn(3, num_points), Xpn(3, num_points);
for( int i = 0; i < num_points; ++i ){
for (int i = 0; i < num_points; ++i) {
Xn.col(i) = X.col(i) - C;
Xpn.col(i) = Xp.col(i) - Cp;
}
// Construct the N matrix (pg. 635).
double Sxx = Xn.row(0).dot(Xpn.row(0));
double Syy = Xn.row(1).dot(Xpn.row(1));
@@ -101,7 +102,7 @@ void AbsoluteOrientation(const Mat3X &X,
Syz - Szy, Sxx - Syy - Szz, Sxy + Syx, Szx + Sxz,
Szx - Sxz, Sxy + Syx, -Sxx + Syy - Szz, Syz + Szy,
Sxy - Syx, Szx + Sxz, Syz + Szy, -Sxx - Syy + Szz;
// Find the unit quaternion q that maximizes qNq. It is the eigenvector
// corresponding to the lagest eigenvalue.
Vec4 q = N.jacobiSvd(Eigen::ComputeFullU).matrixU().col(0);
@@ -180,12 +181,12 @@ static Vec MatrixToConstraint(const Mat &A,
C.setZero();
int idx = 0;
for (int i = 0; i < num_lambda; ++i) {
for( int j = i; j < num_lambda; ++j) {
for (int j = i; j < num_lambda; ++j) {
C(idx) = A(i, j);
if (i != j){
if (i != j) {
C(idx) += A(j, i);
}
++ idx;
++idx;
}
}
return C;
@@ -194,14 +195,14 @@ static Vec MatrixToConstraint(const Mat &A,
// Normalizes the columns of vectors.
static void NormalizeColumnVectors(Mat3X *vectors) {
int num_columns = vectors->cols();
for (int i = 0; i < num_columns; ++i){
for (int i = 0; i < num_columns; ++i) {
vectors->col(i).normalize();
}
}
void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
const Mat3X &X_world,
Mat3 *R,
void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
const Mat3X &X_world,
Mat3 *R,
Vec3 *t) {
CHECK(x_camera.cols() == X_world.cols());
CHECK(x_camera.cols() > 3);
@@ -214,7 +215,7 @@ void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
x_camera_unit.block(0, 0, 2, num_points) = x_camera;
x_camera_unit.row(2).setOnes();
NormalizeColumnVectors(&x_camera_unit);
int num_m_rows = num_points * (num_points - 1) / 2;
int num_tt_variables = num_points * (num_points + 1) / 2;
int num_m_columns = num_tt_variables + 1;
@@ -224,7 +225,7 @@ void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
// Create the constraint equations for the t_ij variables (7) and arrange
// them into the M matrix (8). Also store the initial (i, j) indices.
int row=0;
int row = 0;
for (int i = 0; i < num_points; ++i) {
for (int j = i+1; j < num_points; ++j) {
M(row, row) = -2 * x_camera_unit.col(i).dot(x_camera_unit.col(j));
@@ -243,7 +244,7 @@ void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
}
int num_lambda = num_points + 1; // Dimension of the null space of M.
Mat V = M.jacobiSvd(Eigen::ComputeFullV).matrixV().block(0,
Mat V = M.jacobiSvd(Eigen::ComputeFullV).matrixV().block(0,
num_m_rows,
num_m_columns,
num_lambda);
@@ -263,12 +264,11 @@ void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
int counter_k_row = 0;
for (int idx1 = num_m_rows; idx1 < num_tt_variables; ++idx1) {
for (int idx2 = 0; idx2 < num_m_rows; ++idx2) {
unsigned int i = ij_index(idx1, 0);
unsigned int j = ij_index(idx2, 0);
unsigned int k = ij_index(idx2, 1);
if( i != j && i != k ){
if (i != j && i != k) {
int idx3 = IJToPointIndex(i, j, num_points);
int idx4 = IJToPointIndex(i, k, num_points);
@@ -318,13 +318,12 @@ void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
L_sq = L_sq * Sign(L_sq(IJToIndex(max_L_sq_index,
max_L_sq_index,
num_lambda)));
Vec L(num_lambda);
L(max_L_sq_index) = sqrt(L_sq(IJToIndex(max_L_sq_index,
max_L_sq_index,
num_lambda)));
for (int i = 0; i < num_lambda; ++i) {
if (i != max_L_sq_index) {
L(i) = L_sq(IJToIndex(max_L_sq_index, i, num_lambda)) / L(max_L_sq_index);
@@ -334,7 +333,7 @@ void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
// Correct the scale using the fact that the last constraint is equal to 1.
L = L / (V.row(num_m_columns - 1).dot(L));
Vec X = V * L;
// Recover the distances from the camera center to the 3D points Q.
Vec d(num_points);
d.setZero();
@@ -344,7 +343,7 @@ void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
// Create the 3D points in the camera system.
Mat X_cam(3, num_points);
for (int c_point = 0; c_point < num_points; ++c_point ) {
for (int c_point = 0; c_point < num_points; ++c_point) {
X_cam.col(c_point) = d(c_point) * x_camera_unit.col(c_point);
}
// Recover the camera translation and rotation.
@@ -363,17 +362,17 @@ static void SelectControlPoints(const Mat3X &X_world,
X_control_points->col(0) = mean;
// Computes PCA
X_centered->resize (3, num_points);
X_centered->resize(3, num_points);
for (size_t c = 0; c < num_points; c++) {
X_centered->col(c) = X_world.col (c) - mean;
X_centered->col(c) = X_world.col(c) - mean;
}
Mat3 X_centered_sq = (*X_centered) * X_centered->transpose();
Eigen::JacobiSVD<Mat3> X_centered_sq_svd(X_centered_sq, Eigen::ComputeFullU);
Vec3 w = X_centered_sq_svd.singularValues();
Mat3 u = X_centered_sq_svd.matrixU();
for (size_t c = 0; c < 3; c++) {
double k = sqrt (w (c) / num_points);
X_control_points->col (c + 1) = mean + k * u.col (c);
double k = sqrt(w(c) / num_points);
X_control_points->col(c + 1) = mean + k * u.col(c);
}
}
@@ -382,7 +381,7 @@ static void ComputeBarycentricCoordinates(const Mat3X &X_world_centered,
const Mat34 &X_control_points,
Mat4X *alphas) {
size_t num_points = X_world_centered.cols();
Mat3 C2 ;
Mat3 C2;
for (size_t c = 1; c < 4; c++) {
C2.col(c-1) = X_control_points.col(c) - X_control_points.col(0);
}
@@ -400,7 +399,7 @@ static void ComputeBarycentricCoordinates(const Mat3X &X_world_centered,
// Estimates the coordinates of all real points in the camera coordinate frame
static void ComputePointsCoordinatesInCameraFrame(
const Mat4X &alphas,
const Mat4X &alphas,
const Vec4 &betas,
const Eigen::Matrix<double, 12, 12> &U,
Mat3X *X_camera) {
@@ -423,7 +422,7 @@ static void ComputePointsCoordinatesInCameraFrame(
// Check the sign of the z coordinate of the points (should be positive)
uint num_z_neg = 0;
for (size_t i = 0; i < X_camera->cols(); ++i) {
if ((*X_camera)(2,i) < 0) {
if ((*X_camera)(2, i) < 0) {
num_z_neg++;
}
}
@@ -432,7 +431,7 @@ static void ComputePointsCoordinatesInCameraFrame(
if (num_z_neg > 0.5 * X_camera->cols()) {
C2b = -C2b;
*X_camera = -(*X_camera);
}
}
}
bool EuclideanResectionEPnP(const Mat2X &x_camera,
@@ -442,16 +441,16 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
CHECK(x_camera.cols() == X_world.cols());
CHECK(x_camera.cols() > 3);
size_t num_points = X_world.cols();
// Select the control points.
Mat34 X_control_points;
Mat X_centered;
SelectControlPoints(X_world, &X_centered, &X_control_points);
// Compute the barycentric coordinates.
Mat4X alphas(4, num_points);
ComputeBarycentricCoordinates(X_centered, X_control_points, &alphas);
// Estimates the M matrix with the barycentric coordinates
Mat M(2 * num_points, 12);
Eigen::Matrix<double, 2, 12> sub_M;
@@ -462,17 +461,17 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
double a3 = alphas(3, c);
double ui = x_camera(0, c);
double vi = x_camera(1, c);
M.block(2*c, 0, 2, 12) << a0, 0,
M.block(2*c, 0, 2, 12) << a0, 0,
a0*(-ui), a1, 0,
a1*(-ui), a2, 0,
a1*(-ui), a2, 0,
a2*(-ui), a3, 0,
a3*(-ui), 0,
a3*(-ui), 0,
a0, a0*(-vi), 0,
a1, a1*(-vi), 0,
a2, a2*(-vi), 0,
a3, a3*(-vi);
}
// TODO(julien): Avoid the transpose by rewriting the u2.block() calls.
Eigen::JacobiSVD<Mat> MtMsvd(M.transpose()*M, Eigen::ComputeFullU);
Eigen::Matrix<double, 12, 12> u2 = MtMsvd.matrixU().transpose();
@@ -495,18 +494,18 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
dv2.row(3) = u2.block(10, 3, 1, 3) - u2.block(10, 6, 1, 3);
dv2.row(4) = u2.block(10, 3, 1, 3) - u2.block(10, 9, 1, 3);
dv2.row(5) = u2.block(10, 6, 1, 3) - u2.block(10, 9, 1, 3);
dv3.row(0) = u2.block( 9, 0, 1, 3) - u2.block( 9, 3, 1, 3);
dv3.row(1) = u2.block( 9, 0, 1, 3) - u2.block( 9, 6, 1, 3);
dv3.row(2) = u2.block( 9, 0, 1, 3) - u2.block( 9, 9, 1, 3);
dv3.row(3) = u2.block( 9, 3, 1, 3) - u2.block( 9, 6, 1, 3);
dv3.row(4) = u2.block( 9, 3, 1, 3) - u2.block( 9, 9, 1, 3);
dv3.row(5) = u2.block( 9, 6, 1, 3) - u2.block( 9, 9, 1, 3);
dv4.row(0) = u2.block( 8, 0, 1, 3) - u2.block( 8, 3, 1, 3);
dv4.row(1) = u2.block( 8, 0, 1, 3) - u2.block( 8, 6, 1, 3);
dv4.row(2) = u2.block( 8, 0, 1, 3) - u2.block( 8, 9, 1, 3);
dv4.row(3) = u2.block( 8, 3, 1, 3) - u2.block( 8, 6, 1, 3);
dv4.row(4) = u2.block( 8, 3, 1, 3) - u2.block( 8, 9, 1, 3);
dv4.row(5) = u2.block( 8, 6, 1, 3) - u2.block( 8, 9, 1, 3);
dv3.row(0) = u2.block(9, 0, 1, 3) - u2.block(9, 3, 1, 3);
dv3.row(1) = u2.block(9, 0, 1, 3) - u2.block(9, 6, 1, 3);
dv3.row(2) = u2.block(9, 0, 1, 3) - u2.block(9, 9, 1, 3);
dv3.row(3) = u2.block(9, 3, 1, 3) - u2.block(9, 6, 1, 3);
dv3.row(4) = u2.block(9, 3, 1, 3) - u2.block(9, 9, 1, 3);
dv3.row(5) = u2.block(9, 6, 1, 3) - u2.block(9, 9, 1, 3);
dv4.row(0) = u2.block(8, 0, 1, 3) - u2.block(8, 3, 1, 3);
dv4.row(1) = u2.block(8, 0, 1, 3) - u2.block(8, 6, 1, 3);
dv4.row(2) = u2.block(8, 0, 1, 3) - u2.block(8, 9, 1, 3);
dv4.row(3) = u2.block(8, 3, 1, 3) - u2.block(8, 6, 1, 3);
dv4.row(4) = u2.block(8, 3, 1, 3) - u2.block(8, 9, 1, 3);
dv4.row(5) = u2.block(8, 6, 1, 3) - u2.block(8, 9, 1, 3);
Eigen::Matrix<double, 6, 10> L;
for (size_t r = 0; r < 6; r++) {
@@ -520,7 +519,7 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
2.0 * dv2.row(r).dot(dv4.row(r)),
2.0 * dv3.row(r).dot(dv4.row(r)),
dv4.row(r).dot(dv4.row(r));
}
}
Vec6 rho;
rho << (X_control_points.col(0) - X_control_points.col(1)).squaredNorm(),
(X_control_points.col(0) - X_control_points.col(2)).squaredNorm(),
@@ -528,7 +527,7 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
(X_control_points.col(1) - X_control_points.col(2)).squaredNorm(),
(X_control_points.col(1) - X_control_points.col(3)).squaredNorm(),
(X_control_points.col(2) - X_control_points.col(3)).squaredNorm();
// There are three possible solutions based on the three approximations of L
// (betas). Below, each one is solved for then the best one is chosen.
Mat3X X_camera;
@@ -548,7 +547,8 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
// check, is the right approach. So far this seems the case.
//
// TODO(sergey): Made it an option for now, in some cases it makes sense to
// still fallback to reprojection solution (see bug [#32765] from Blender bug tracker)
// still fallback to reprojection solution
// For details see bug [#32765] from Blender bug tracker
// double kSuccessThreshold = std::numeric_limits<double>::max();
double kSuccessThreshold = success_threshold;
@@ -559,15 +559,15 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
Vec4 betas = Vec4::Zero();
Eigen::Matrix<double, 6, 4> l_6x4;
for (size_t r = 0; r < 6; r++) {
l_6x4.row(r) << L(r, 0), L(r, 1), L(r, 3), L(r, 6);
l_6x4.row(r) << L(r, 0), L(r, 1), L(r, 3), L(r, 6);
}
Eigen::JacobiSVD<Mat> svd_of_l4(l_6x4,
Eigen::JacobiSVD<Mat> svd_of_l4(l_6x4,
Eigen::ComputeFullU | Eigen::ComputeFullV);
Vec4 b4 = svd_of_l4.solve(rho);
if ((l_6x4 * b4).isApprox(rho, kSuccessThreshold)) {
if (b4(0) < 0) {
b4 = -b4;
}
}
b4(0) = std::sqrt(b4(0));
betas << b4(0), b4(1) / b4(0), b4(2) / b4(0), b4(3) / b4(0);
ComputePointsCoordinatesInCameraFrame(alphas, betas, u2, &X_camera);
@@ -578,14 +578,14 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
ts[0].setZero();
rmse(0) = std::numeric_limits<double>::max();
}
// Find the second possible solution for R, t corresponding to:
// Betas = [b00 b01 b11 b02 b12 b22 b03 b13 b23 b33]
// Betas_approx_2 = [b00 b01 b11]
betas.setZero();
Eigen::Matrix<double, 6, 3> l_6x3;
l_6x3 = L.block(0, 0, 6, 3);
Eigen::JacobiSVD<Mat> svdOfL3(l_6x3,
Eigen::JacobiSVD<Mat> svdOfL3(l_6x3,
Eigen::ComputeFullU | Eigen::ComputeFullV);
Vec3 b3 = svdOfL3.solve(rho);
VLOG(2) << " rho = " << rho;
@@ -611,14 +611,14 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
ts[1].setZero();
rmse(1) = std::numeric_limits<double>::max();
}
// Find the third possible solution for R, t corresponding to:
// Betas = [b00 b01 b11 b02 b12 b22 b03 b13 b23 b33]
// Betas_approx_3 = [b00 b01 b11 b02 b12]
betas.setZero();
Eigen::Matrix<double, 6, 5> l_6x5;
l_6x5 = L.block(0, 0, 6, 5);
Eigen::JacobiSVD<Mat> svdOfL5(l_6x5,
Eigen::JacobiSVD<Mat> svdOfL5(l_6x5,
Eigen::ComputeFullU | Eigen::ComputeFullV);
Vec5 b5 = svdOfL5.solve(rho);
if ((l_6x5 * b5).isApprox(rho, kSuccessThreshold)) {
@@ -650,7 +650,7 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
ts[2].setZero();
rmse(2) = std::numeric_limits<double>::max();
}
// Finally, with all three solutions, select the (R, t) with the best RMSE.
VLOG(2) << "RMSE for solution 0: " << rmse(0);
VLOG(2) << "RMSE for solution 1: " << rmse(1);
@@ -675,5 +675,5 @@ bool EuclideanResectionEPnP(const Mat2X &x_camera,
return true;
}
} // namespace resection
} // namespace libmv
} // namespace resection
} // namespace libmv

14
extern/libmv/libmv/multiview/euclidean_resection.h vendored
View File

@@ -26,7 +26,7 @@
namespace libmv {
namespace euclidean_resection {
enum ResectionMethod {
RESECTION_ANSAR_DANIILIDIS,
@@ -48,7 +48,7 @@ enum ResectionMethod {
* \param success_threshold Threshold of an error which is still considered a success
* (currently used by EPnP algorithm only)
*/
bool EuclideanResection(const Mat2X &x_camera,
bool EuclideanResection(const Mat2X &x_camera,
const Mat3X &X_world,
Mat3 *R, Vec3 *t,
ResectionMethod method = RESECTION_EPNP,
@@ -67,7 +67,7 @@ bool EuclideanResection(const Mat2X &x_camera,
* \param t Solution for the camera translation vector
* \param method Resection method
*/
bool EuclideanResection(const Mat &x_image,
bool EuclideanResection(const Mat &x_image,
const Mat3X &X_world,
const Mat3 &K,
Mat3 *R, Vec3 *t,
@@ -101,7 +101,7 @@ void AbsoluteOrientation(const Mat3X &X,
* This is the algorithm described in: "Linear Pose Estimation from Points or
* Lines", by Ansar, A. and Daniilidis, PAMI 2003. vol. 25, no. 5.
*/
void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
const Mat3X &X_world,
Mat3 *R, Vec3 *t);
/**
@@ -121,12 +121,12 @@ void EuclideanResectionAnsarDaniilidis(const Mat2X &x_camera,
* \note: the non-linear optimization is not implemented here.
*/
bool EuclideanResectionEPnP(const Mat2X &x_camera,
const Mat3X &X_world,
const Mat3X &X_world,
Mat3 *R, Vec3 *t,
double success_threshold = 1e-3);
} // namespace euclidean_resection
} // namespace libmv
} // namespace euclidean_resection
} // namespace libmv
#endif /* LIBMV_MULTIVIEW_EUCLIDEAN_RESECTION_H_ */

4
extern/libmv/libmv/multiview/fundamental.cc vendored
View File

@@ -18,13 +18,14 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/multiview/fundamental.h"
#include "libmv/logging/logging.h"
#include "libmv/numeric/numeric.h"
#include "libmv/numeric/poly.h"
#include "libmv/multiview/conditioning.h"
#include "libmv/multiview/projection.h"
#include "libmv/multiview/triangulation.h"
#include "libmv/multiview/fundamental.h"
namespace libmv {
@@ -304,7 +305,6 @@ void MotionFromEssential(const Mat3 &E,
std::vector<Vec3> *ts) {
Eigen::JacobiSVD<Mat3> USV(E, Eigen::ComputeFullU | Eigen::ComputeFullV);
Mat3 U = USV.matrixU();
Vec3 d = USV.singularValues();
Mat3 Vt = USV.matrixV().transpose();
// Last column of U is undetermined since d = (a a 0).

189
extern/libmv/libmv/multiview/homography.cc vendored
View File

@@ -18,8 +18,9 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/logging/logging.h"
#include "libmv/multiview/homography.h"
#include "libmv/logging/logging.h"
#include "libmv/multiview/homography_parameterization.h"
namespace libmv {
@@ -55,30 +56,30 @@ static bool Homography2DFromCorrespondencesLinearEuc(
Mat b = Mat::Zero(n * 3, 1);
for (int i = 0; i < n; ++i) {
int j = 3 * i;
L(j, 0) = x1(0, i); // a
L(j, 1) = x1(1, i); // b
L(j, 2) = 1.0; // c
L(j, 6) = -x2(0, i) * x1(0, i); // g
L(j, 7) = -x2(0, i) * x1(1, i); // h
b(j, 0) = x2(0, i); // i
L(j, 0) = x1(0, i); // a
L(j, 1) = x1(1, i); // b
L(j, 2) = 1.0; // c
L(j, 6) = -x2(0, i) * x1(0, i); // g
L(j, 7) = -x2(0, i) * x1(1, i); // h
b(j, 0) = x2(0, i); // i
++j;
L(j, 3) = x1(0, i); // d
L(j, 4) = x1(1, i); // e
L(j, 5) = 1.0; // f
L(j, 6) = -x2(1, i) * x1(0, i); // g
L(j, 7) = -x2(1, i) * x1(1, i); // h
b(j, 0) = x2(1, i); // i
L(j, 3) = x1(0, i); // d
L(j, 4) = x1(1, i); // e
L(j, 5) = 1.0; // f
L(j, 6) = -x2(1, i) * x1(0, i); // g
L(j, 7) = -x2(1, i) * x1(1, i); // h
b(j, 0) = x2(1, i); // i
// This ensures better stability
// TODO(julien) make a lite version without this 3rd set
++j;
L(j, 0) = x2(1, i) * x1(0, i); // a
L(j, 1) = x2(1, i) * x1(1, i); // b
L(j, 2) = x2(1, i); // c
L(j, 3) = -x2(0, i) * x1(0, i); // d
L(j, 4) = -x2(0, i) * x1(1, i); // e
L(j, 5) = -x2(0, i) ; // f
L(j, 0) = x2(1, i) * x1(0, i); // a
L(j, 1) = x2(1, i) * x1(1, i); // b
L(j, 2) = x2(1, i); // c
L(j, 3) = -x2(0, i) * x1(0, i); // d
L(j, 4) = -x2(0, i) * x1(1, i); // e
L(j, 5) = -x2(0, i); // f
}
// Solve Lx=B
Vec h = L.fullPivLu().solve(b);
@@ -103,7 +104,7 @@ bool Homography2DFromCorrespondencesLinear(const Mat &x1,
Mat3 *H,
double expected_precision) {
if (x1.rows() == 2) {
return Homography2DFromCorrespondencesLinearEuc(x1, x2, H,
return Homography2DFromCorrespondencesLinearEuc(x1, x2, H,
expected_precision);
}
assert(3 == x1.rows());
@@ -119,29 +120,29 @@ bool Homography2DFromCorrespondencesLinear(const Mat &x1,
Mat b = Mat::Zero(n * 3, 1);
for (int i = 0; i < n; ++i) {
int j = 3 * i;
L(j, 0) = x2(w, i) * x1(x, i);//a
L(j, 1) = x2(w, i) * x1(y, i);//b
L(j, 2) = x2(w, i) * x1(w, i);//c
L(j, 6) = -x2(x, i) * x1(x, i);//g
L(j, 7) = -x2(x, i) * x1(y, i);//h
L(j, 0) = x2(w, i) * x1(x, i); // a
L(j, 1) = x2(w, i) * x1(y, i); // b
L(j, 2) = x2(w, i) * x1(w, i); // c
L(j, 6) = -x2(x, i) * x1(x, i); // g
L(j, 7) = -x2(x, i) * x1(y, i); // h
b(j, 0) = x2(x, i) * x1(w, i);
++j;
L(j, 3) = x2(w, i) * x1(x, i);//d
L(j, 4) = x2(w, i) * x1(y, i);//e
L(j, 5) = x2(w, i) * x1(w, i);//f
L(j, 6) = -x2(y, i) * x1(x, i);//g
L(j, 7) = -x2(y, i) * x1(y, i);//h
L(j, 3) = x2(w, i) * x1(x, i); // d
L(j, 4) = x2(w, i) * x1(y, i); // e
L(j, 5) = x2(w, i) * x1(w, i); // f
L(j, 6) = -x2(y, i) * x1(x, i); // g
L(j, 7) = -x2(y, i) * x1(y, i); // h
b(j, 0) = x2(y, i) * x1(w, i);
// This ensures better stability
++j;
L(j, 0) = x2(y, i) * x1(x, i);//a
L(j, 1) = x2(y, i) * x1(y, i);//b
L(j, 2) = x2(y, i) * x1(w, i);//c
L(j, 3) = -x2(x, i) * x1(x, i);//d
L(j, 4) = -x2(x, i) * x1(y, i);//e
L(j, 5) = -x2(x, i) * x1(w, i);//f
L(j, 0) = x2(y, i) * x1(x, i); // a
L(j, 1) = x2(y, i) * x1(y, i); // b
L(j, 2) = x2(y, i) * x1(w, i); // c
L(j, 3) = -x2(x, i) * x1(x, i); // d
L(j, 4) = -x2(x, i) * x1(y, i); // e
L(j, 5) = -x2(x, i) * x1(w, i); // f
}
// Solve Lx=B
Vec h = L.fullPivLu().solve(b);
@@ -165,11 +166,11 @@ bool Homography2DFromCorrespondencesLinear(const Mat &x1,
* |0 0 0 1| -> |-x2w| |-1 0 0 0| -> | x2x| |0 0 0 0| -> | 0 |
* |0 0 0 0| | 0 | | 0 0 0 0| | 0 | |0 0 0 1| |-x2w|
* |0-1 0 0| | x2y| | 0 0 0 0| | 0 | |0 0-1 0| | x2z|
* |a b c d|
* A = |e f g h|
* |a b c d|
* A = |e f g h|
* |i j k l|
* |m n o 1|
*
* |m n o 1|
*
* x2^t * H1 * A *x1 = (-x2w*a +x2x*m )*x1x + (-x2w*b +x2x*n )*x1y + (-x2w*c +x2x*o )*x1z + (-x2w*d +x2x*1 )*x1w = 0
* x2^t * H2 * A *x1 = (-x2z*e +x2y*i )*x1x + (-x2z*f +x2y*j )*x1y + (-x2z*g +x2y*k )*x1z + (-x2z*h +x2y*l )*x1w = 0
* x2^t * H3 * A *x1 = (-x2z*a +x2x*i )*x1x + (-x2z*b +x2x*j )*x1y + (-x2z*c +x2x*k )*x1z + (-x2z*d +x2x*l )*x1w = 0
@@ -196,64 +197,64 @@ bool Homography3DFromCorrespondencesLinear(const Mat &x1,
Mat b = Mat::Zero(n * 6, 1);
for (int i = 0; i < n; ++i) {
int j = 6 * i;
L(j, 0) = -x2(w, i) * x1(x, i);//a
L(j, 1) = -x2(w, i) * x1(y, i);//b
L(j, 2) = -x2(w, i) * x1(z, i);//c
L(j, 3) = -x2(w, i) * x1(w, i);//d
L(j,12) = x2(x, i) * x1(x, i);//m
L(j,13) = x2(x, i) * x1(y, i);//n
L(j,14) = x2(x, i) * x1(z, i);//o
b(j, 0) = -x2(x, i) * x1(w, i);
L(j, 0) = -x2(w, i) * x1(x, i); // a
L(j, 1) = -x2(w, i) * x1(y, i); // b
L(j, 2) = -x2(w, i) * x1(z, i); // c
L(j, 3) = -x2(w, i) * x1(w, i); // d
L(j, 12) = x2(x, i) * x1(x, i); // m
L(j, 13) = x2(x, i) * x1(y, i); // n
L(j, 14) = x2(x, i) * x1(z, i); // o
b(j, 0) = -x2(x, i) * x1(w, i);
++j;
L(j, 4) = -x2(z, i) * x1(x, i);//e
L(j, 5) = -x2(z, i) * x1(y, i);//f
L(j, 6) = -x2(z, i) * x1(z, i);//g
L(j, 7) = -x2(z, i) * x1(w, i);//h
L(j, 8) = x2(y, i) * x1(x, i);//i
L(j, 9) = x2(y, i) * x1(y, i);//j
L(j,10) = x2(y, i) * x1(z, i);//k
L(j,11) = x2(y, i) * x1(w, i);//l
L(j, 4) = -x2(z, i) * x1(x, i); // e
L(j, 5) = -x2(z, i) * x1(y, i); // f
L(j, 6) = -x2(z, i) * x1(z, i); // g
L(j, 7) = -x2(z, i) * x1(w, i); // h
L(j, 8) = x2(y, i) * x1(x, i); // i
L(j, 9) = x2(y, i) * x1(y, i); // j
L(j, 10) = x2(y, i) * x1(z, i); // k
L(j, 11) = x2(y, i) * x1(w, i); // l
++j;
L(j, 0) = -x2(z, i) * x1(x, i);//a
L(j, 1) = -x2(z, i) * x1(y, i);//b
L(j, 2) = -x2(z, i) * x1(z, i);//c
L(j, 3) = -x2(z, i) * x1(w, i);//d
L(j, 8) = x2(x, i) * x1(x, i);//i
L(j, 9) = x2(x, i) * x1(y, i);//j
L(j,10) = x2(x, i) * x1(z, i);//k
L(j,11) = x2(x, i) * x1(w, i);//l
L(j, 0) = -x2(z, i) * x1(x, i); // a
L(j, 1) = -x2(z, i) * x1(y, i); // b
L(j, 2) = -x2(z, i) * x1(z, i); // c
L(j, 3) = -x2(z, i) * x1(w, i); // d
L(j, 8) = x2(x, i) * x1(x, i); // i
L(j, 9) = x2(x, i) * x1(y, i); // j
L(j, 10) = x2(x, i) * x1(z, i); // k
L(j, 11) = x2(x, i) * x1(w, i); // l
++j;
L(j, 4) = -x2(w, i) * x1(x, i);//e
L(j, 5) = -x2(w, i) * x1(y, i);//f
L(j, 6) = -x2(w, i) * x1(z, i);//g
L(j, 7) = -x2(w, i) * x1(w, i);//h
L(j,12) = x2(y, i) * x1(x, i);//m
L(j,13) = x2(y, i) * x1(y, i);//n
L(j,14) = x2(y, i) * x1(z, i);//o
b(j, 0) = -x2(y, i) * x1(w, i);
L(j, 4) = -x2(w, i) * x1(x, i); // e
L(j, 5) = -x2(w, i) * x1(y, i); // f
L(j, 6) = -x2(w, i) * x1(z, i); // g
L(j, 7) = -x2(w, i) * x1(w, i); // h
L(j, 12) = x2(y, i) * x1(x, i); // m
L(j, 13) = x2(y, i) * x1(y, i); // n
L(j, 14) = x2(y, i) * x1(z, i); // o
b(j, 0) = -x2(y, i) * x1(w, i);
++j;
L(j, 0) = -x2(y, i) * x1(x, i);//a
L(j, 1) = -x2(y, i) * x1(y, i);//b
L(j, 2) = -x2(y, i) * x1(z, i);//c
L(j, 3) = -x2(y, i) * x1(w, i);//d
L(j, 4) = x2(x, i) * x1(x, i);//e
L(j, 5) = x2(x, i) * x1(y, i);//f
L(j, 6) = x2(x, i) * x1(z, i);//g
L(j, 7) = x2(x, i) * x1(w, i);//h
L(j, 0) = -x2(y, i) * x1(x, i); // a
L(j, 1) = -x2(y, i) * x1(y, i); // b
L(j, 2) = -x2(y, i) * x1(z, i); // c
L(j, 3) = -x2(y, i) * x1(w, i); // d
L(j, 4) = x2(x, i) * x1(x, i); // e
L(j, 5) = x2(x, i) * x1(y, i); // f
L(j, 6) = x2(x, i) * x1(z, i); // g
L(j, 7) = x2(x, i) * x1(w, i); // h
++j;
L(j, 8) = -x2(w, i) * x1(x, i);//i
L(j, 9) = -x2(w, i) * x1(y, i);//j
L(j,10) = -x2(w, i) * x1(z, i);//k
L(j,11) = -x2(w, i) * x1(w, i);//l
L(j,12) = x2(z, i) * x1(x, i);//m
L(j,13) = x2(z, i) * x1(y, i);//n
L(j,14) = x2(z, i) * x1(z, i);//o
b(j, 0) = -x2(z, i) * x1(w, i);
L(j, 8) = -x2(w, i) * x1(x, i); // i
L(j, 9) = -x2(w, i) * x1(y, i); // j
L(j, 10) = -x2(w, i) * x1(z, i); // k
L(j, 11) = -x2(w, i) * x1(w, i); // l
L(j, 12) = x2(z, i) * x1(x, i); // m
L(j, 13) = x2(z, i) * x1(y, i); // n
L(j, 14) = x2(z, i) * x1(z, i); // o
b(j, 0) = -x2(z, i) * x1(w, i);
}
// Solve Lx=B
Vec h = L.fullPivLu().solve(b);

12
extern/libmv/libmv/multiview/homography.h vendored
View File

@@ -50,12 +50,12 @@ namespace libmv {
bool Homography2DFromCorrespondencesLinear(const Mat &x1,
const Mat &x2,
Mat3 *H,
double expected_precision =
double expected_precision =
EigenDouble::dummy_precision());
/**
* 3D Homography transformation estimation.
*
*
* This function can be used in order to estimate the homography transformation
* from a list of 3D correspondences.
*
@@ -68,16 +68,16 @@ bool Homography2DFromCorrespondencesLinear(const Mat &x1,
* |m n o 1|
* \param[in] expected_precision The expected precision in order for instance
* to accept almost homography matrices.
*
*
* \return true if the transformation estimation has succeeded
*
*
* \note Need at least 5 non coplanar points
* \note Points coordinates must be in homogeneous coordinates
*/
bool Homography3DFromCorrespondencesLinear(const Mat &x1,
const Mat &x2,
Mat4 *H,
double expected_precision =
double expected_precision =
EigenDouble::dummy_precision());
/**
@@ -87,6 +87,6 @@ bool Homography3DFromCorrespondencesLinear(const Mat &x1,
*/
double SymmetricGeometricDistance(Mat3 &H, Vec2 &x1, Vec2 &x2);
} // namespace libmv
} // namespace libmv
#endif // LIBMV_MULTIVIEW_HOMOGRAPHY_H_

12
extern/libmv/libmv/multiview/homography_parameterization.h vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2011 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -40,7 +40,7 @@ class Homography2DNormalizedParameterization {
typedef Eigen::Matrix<T, 3, 3> Parameterized; // H
/// Convert from the 8 parameters to a H matrix.
static void To(const Parameters &p, Parameterized *h) {
static void To(const Parameters &p, Parameterized *h) {
*h << p(0), p(1), p(2),
p(3), p(4), p(5),
p(6), p(7), 1.0;
@@ -70,7 +70,7 @@ class Homography3DNormalizedParameterization {
typedef Eigen::Matrix<T, 4, 4> Parameterized; // H
/// Convert from the 15 parameters to a H matrix.
static void To(const Parameters &p, Parameterized *h) {
static void To(const Parameters &p, Parameterized *h) {
*h << p(0), p(1), p(2), p(3),
p(4), p(5), p(6), p(7),
p(8), p(9), p(10), p(11),
@@ -86,6 +86,6 @@ class Homography3DNormalizedParameterization {
}
};
} // namespace libmv
} // namespace libmv
#endif // LIBMV_MULTIVIEW_HOMOGRAPHY_PARAMETERIZATION_H_

8
extern/libmv/libmv/multiview/nviewtriangulation.h vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2009 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -17,7 +17,7 @@
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
//
//
// Compute a 3D position of a point from several images of it. In particular,
// compute the projective point X in R^4 such that x = PX.
//

125
extern/libmv/libmv/multiview/panography.cc vendored Normal file
View File

@@ -0,0 +1,125 @@
// Copyright (c) 2009 libmv authors.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
//
#include "libmv/multiview/panography.h"
namespace libmv {
static bool Build_Minimal2Point_PolynomialFactor(
const Mat & x1, const Mat & x2,
double * P) { // P must be a double[4]
assert(2 == x1.rows());
assert(2 == x1.cols());
assert(x1.rows() == x2.rows());
assert(x1.cols() == x2.cols());
// Setup the variable of the input problem:
Vec xx1 = (x1.col(0)).transpose();
Vec yx1 = (x1.col(1)).transpose();
double a12 = xx1.dot(yx1);
Vec xx2 = (x2.col(0)).transpose();
Vec yx2 = (x2.col(1)).transpose();
double b12 = xx2.dot(yx2);
double a1 = xx1.squaredNorm();
double a2 = yx1.squaredNorm();
double b1 = xx2.squaredNorm();
double b2 = yx2.squaredNorm();
// Build the 3rd degre polynomial in F^2.
//
// f^6 * p + f^4 * q + f^2* r + s = 0;
//
// Coefficients in ascending powers of alpha, i.e. P[N]*x^N.
// Run panography_coeffs.py to get the below coefficients.
P[0] = b1*b2*a12*a12-a1*a2*b12*b12;
P[1] = -2*a1*a2*b12+2*a12*b1*b2+b1*a12*a12+b2*a12*a12-a1*b12*b12-a2*b12*b12;
P[2] = b1*b2-a1*a2-2*a1*b12-2*a2*b12+2*a12*b1+2*a12*b2+a12*a12-b12*b12;
P[3] = b1+b2-2*b12-a1-a2+2*a12;
// If P[3] equal to 0 we get ill conditionned data
return (P[3] != 0.0);
}
// This implements a minimal solution (2 points) for panoramic stitching:
//
// http://www.cs.ubc.ca/~mbrown/minimal/minimal.html
//
// [1] M. Brown and R. Hartley and D. Nister. Minimal Solutions for Panoramic
// Stitching. CVPR07.
void F_FromCorrespondance_2points(const Mat &x1, const Mat &x2,
vector<double> *fs) {
// Build Polynomial factor to get squared focal value.
double P[4];
Build_Minimal2Point_PolynomialFactor(x1, x2, &P[0]);
// Solve it by using F = f^2 and a Cubic polynomial solver
//
// F^3 * p + F^2 * q + F^1 * r + s = 0
//
double roots[3];
int num_roots = SolveCubicPolynomial(P, roots);
for (int i = 0; i < num_roots; ++i) {
if (roots[i] > 0.0) {
fs->push_back(sqrt(roots[i]));
}
}
}
// Compute the 3x3 rotation matrix that fits two 3D point clouds in the least
// square sense. The method is from:
//
// K. Arun,T. Huand and D. Blostein. Least-squares fitting of 2 3-D point
// sets. IEEE Transactions on Pattern Analysis and Machine Intelligence,
// 9:698-700, 1987.
void GetR_FixedCameraCenter(const Mat &x1, const Mat &x2,
const double focal,
Mat3 *R) {
assert(3 == x1.rows());
assert(2 <= x1.cols());
assert(x1.rows() == x2.rows());
assert(x1.cols() == x2.cols());
// Build simplified K matrix
Mat3 K(Mat3::Identity() * 1.0/focal);
K(2, 2)= 1.0;
// Build the correlation matrix; equation (22) in [1].
Mat3 C = Mat3::Zero();
for (int i = 0; i < x1.cols(); ++i) {
Mat r1i = (K * x1.col(i)).normalized();
Mat r2i = (K * x2.col(i)).normalized();
C += r2i * r1i.transpose();
}
// Solve for rotation. Equations (24) and (25) in [1].
Eigen::JacobiSVD<Mat> svd(C, Eigen::ComputeThinU | Eigen::ComputeThinV);
Mat3 scale = Mat3::Identity();
scale(2, 2) = ((svd.matrixU() * svd.matrixV().transpose()).determinant() > 0.0)
? 1.0
: -1.0;
(*R) = svd.matrixU() * scale * svd.matrixV().transpose();
}
} // namespace libmv

94
extern/libmv/libmv/multiview/panography.h vendored
View File

@@ -28,45 +28,6 @@
namespace libmv {
static bool Build_Minimal2Point_PolynomialFactor(
const Mat & x1, const Mat & x2,
double * P) //P must be a double[4]
{
assert(2 == x1.rows());
assert(2 == x1.cols());
assert(x1.rows() == x2.rows());
assert(x1.cols() == x2.cols());
// Setup the variable of the input problem:
Vec xx1 = (x1.col(0)).transpose();
Vec yx1 = (x1.col(1)).transpose();
double a12 = xx1.dot(yx1);
Vec xx2 = (x2.col(0)).transpose();
Vec yx2 = (x2.col(1)).transpose();
double b12 = xx2.dot(yx2);
double a1 = xx1.squaredNorm();
double a2 = yx1.squaredNorm();
double b1 = xx2.squaredNorm();
double b2 = yx2.squaredNorm();
// Build the 3rd degre polynomial in F^2.
//
// f^6 * p + f^4 * q + f^2* r + s = 0;
//
// Coefficients in ascending powers of alpha, i.e. P[N]*x^N.
// Run panography_coeffs.py to get the below coefficients.
P[0] = b1*b2*a12*a12-a1*a2*b12*b12;
P[1] = -2*a1*a2*b12+2*a12*b1*b2+b1*a12*a12+b2*a12*a12-a1*b12*b12-a2*b12*b12;
P[2] = b1*b2-a1*a2-2*a1*b12-2*a2*b12+2*a12*b1+2*a12*b2+a12*a12-b12*b12;
P[3] = b1+b2-2*b12-a1-a2+2*a12;
// If P[3] equal to 0 we get ill conditionned data
return (P[3] != 0.0);
}
// This implements a minimal solution (2 points) for panoramic stitching:
//
// http://www.cs.ubc.ca/~mbrown/minimal/minimal.html
@@ -92,25 +53,8 @@ static bool Build_Minimal2Point_PolynomialFactor(
// K = [0 f 0]
// [0 0 1]
//
static void F_FromCorrespondance_2points(const Mat &x1, const Mat &x2,
vector<double> *fs) {
// Build Polynomial factor to get squared focal value.
double P[4];
Build_Minimal2Point_PolynomialFactor(x1, x2, &P[0]);
// Solve it by using F = f^2 and a Cubic polynomial solver
//
// F^3 * p + F^2 * q + F^1 * r + s = 0
//
double roots[3];
int num_roots = SolveCubicPolynomial(P, roots);
for (int i = 0; i < num_roots; ++i) {
if (roots[i] > 0.0) {
fs->push_back(sqrt(roots[i]));
}
}
}
void F_FromCorrespondance_2points(const Mat &x1, const Mat &x2,
vector<double> *fs);
// Compute the 3x3 rotation matrix that fits two 3D point clouds in the least
// square sense. The method is from:
@@ -146,36 +90,10 @@ static void F_FromCorrespondance_2points(const Mat &x1, const Mat &x2,
//
// R = arg min || X2 - R * x1 ||
//
static void GetR_FixedCameraCenter(const Mat &x1, const Mat &x2,
const double focal,
Mat3 *R) {
assert(3 == x1.rows());
assert(2 <= x1.cols());
assert(x1.rows() == x2.rows());
assert(x1.cols() == x2.cols());
// Build simplified K matrix
Mat3 K( Mat3::Identity() * 1.0/focal );
K(2,2)= 1.0;
// Build the correlation matrix; equation (22) in [1].
Mat3 C = Mat3::Zero();
for(int i = 0; i < x1.cols(); ++i) {
Mat r1i = (K * x1.col(i)).normalized();
Mat r2i = (K * x2.col(i)).normalized();
C += r2i * r1i.transpose();
}
// Solve for rotation. Equations (24) and (25) in [1].
Eigen::JacobiSVD<Mat> svd(C, Eigen::ComputeThinU | Eigen::ComputeThinV);
Mat3 scale = Mat3::Identity();
scale(2,2) = ((svd.matrixU() * svd.matrixV().transpose()).determinant() > 0.0)
? 1.0
: -1.0;
(*R) = svd.matrixU() * scale * svd.matrixV().transpose();
}
void GetR_FixedCameraCenter(const Mat &x1, const Mat &x2,
const double focal,
Mat3 *R);
} // namespace libmv
#endif // LIBMV_MULTIVIEW_PANOGRAPHY_H
#endif // LIBMV_MULTIVIEW_PANOGRAPHY_H

59
extern/libmv/libmv/multiview/projection.cc vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2007, 2008 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -38,10 +38,11 @@ void KRt_From_P(const Mat34 &P, Mat3 *Kp, Mat3 *Rp, Vec3 *tp) {
// Set K(2,1) to zero.
if (K(2, 1) != 0) {
double c = -K(2,2);
double s = K(2,1);
double c = -K(2, 2);
double s = K(2, 1);
double l = sqrt(c * c + s * s);
c /= l; s /= l;
c /= l;
s /= l;
Mat3 Qx;
Qx << 1, 0, 0,
0, c, -s,
@@ -54,7 +55,8 @@ void KRt_From_P(const Mat34 &P, Mat3 *Kp, Mat3 *Rp, Vec3 *tp) {
double c = K(2, 2);
double s = K(2, 0);
double l = sqrt(c * c + s * s);
c /= l; s /= l;
c /= l;
s /= l;
Mat3 Qy;
Qy << c, 0, s,
0, 1, 0,
@@ -67,11 +69,12 @@ void KRt_From_P(const Mat34 &P, Mat3 *Kp, Mat3 *Rp, Vec3 *tp) {
double c = -K(1, 1);
double s = K(1, 0);
double l = sqrt(c * c + s * s);
c /= l; s /= l;
c /= l;
s /= l;
Mat3 Qz;
Qz << c,-s, 0,
s, c, 0,
0, 0, 1;
Qz << c, -s, 0,
s, c, 0,
0, 0, 1;
K = K * Qz;
Q = Qz.transpose() * Q;
}
@@ -83,19 +86,19 @@ void KRt_From_P(const Mat34 &P, Mat3 *Kp, Mat3 *Rp, Vec3 *tp) {
// - K(0,0) > 0
// - K(2,2) = 1
// - det(R) = 1
if (K(2,2) < 0) {
if (K(2, 2) < 0) {
K = -K;
R = -R;
}
if (K(1,1) < 0) {
if (K(1, 1) < 0) {
Mat3 S;
S << 1, 0, 0,
0,-1, 0,
0, 0, 1;
S << 1, 0, 0,
0, -1, 0,
0, 0, 1;
K = K * S;
R = S * R;
}
if (K(0,0) < 0) {
if (K(0, 0) < 0) {
Mat3 S;
S << -1, 0, 0,
0, 1, 0,
@@ -106,13 +109,13 @@ void KRt_From_P(const Mat34 &P, Mat3 *Kp, Mat3 *Rp, Vec3 *tp) {
// Compute translation.
Vec p(3);
p << P(0,3), P(1,3), P(2,3);
p << P(0, 3), P(1, 3), P(2, 3);
// TODO(pau) This sould be done by a SolveLinearSystem(A, b, &x) call.
// TODO(keir) use the eigen LU solver syntax...
Vec3 t = K.inverse() * p;
// scale K so that K(2,2) = 1
K = K / K(2,2);
K = K / K(2, 2);
*Kp = K;
*Rp = R;
@@ -140,10 +143,10 @@ void ProjectionChangeAspectRatio(const Mat34 &P,
0, aspect_ratio_new / aspect_ratio, 0,
0, 0, 1;
Mat34 P_temp;
ProjectionShiftPrincipalPoint(P, principal_point, Vec2(0,0), &P_temp);
ProjectionShiftPrincipalPoint(P, principal_point, Vec2(0, 0), &P_temp);
P_temp = T * P_temp;
ProjectionShiftPrincipalPoint(P_temp, Vec2(0,0), principal_point, P_new);
ProjectionShiftPrincipalPoint(P_temp, Vec2(0, 0), principal_point, P_new);
}
void HomogeneousToEuclidean(const Mat &H, Mat *X) {
@@ -198,15 +201,15 @@ void EuclideanToHomogeneous(const Vec3 &X, Vec4 *H) {
// TODO(julien) Call conditioning.h/ApplyTransformationToPoints ?
void EuclideanToNormalizedCamera(const Mat2X &x, const Mat3 &K, Mat2X *n) {
Mat3X x_image_h;
EuclideanToHomogeneous(x, &x_image_h);
Mat3X x_camera_h = K.inverse() * x_image_h;
HomogeneousToEuclidean(x_camera_h, n);
Mat3X x_image_h;
EuclideanToHomogeneous(x, &x_image_h);
Mat3X x_camera_h = K.inverse() * x_image_h;
HomogeneousToEuclidean(x_camera_h, n);
}
void HomogeneousToNormalizedCamera(const Mat3X &x, const Mat3 &K, Mat2X *n) {
Mat3X x_camera_h = K.inverse() * x;
HomogeneousToEuclidean(x_camera_h, n);
Mat3X x_camera_h = K.inverse() * x;
HomogeneousToEuclidean(x_camera_h, n);
}
double Depth(const Mat3 &R, const Vec3 &t, const Vec3 &X) {

24
extern/libmv/libmv/multiview/projection.h vendored
View File

@@ -94,7 +94,7 @@ inline Vec3 EuclideanToHomogeneous(const Vec2 &x) {
inline Vec4 EuclideanToHomogeneous(const Vec3 &x) {
return Vec4(x(0), x(1), x(2), 1);
}
// Conversion from image coordinates to normalized camera coordinates
// Conversion from image coordinates to normalized camera coordinates
void EuclideanToNormalizedCamera(const Mat2X &x, const Mat3 &K, Mat2X *n);
void HomogeneousToNormalizedCamera(const Mat3X &x, const Mat3 &K, Mat2X *n);
@@ -180,16 +180,16 @@ double Depth(const Mat3 &R, const Vec3 &t, const Vec4 &X);
* Returns true if the homogenious 3D point X is in front of
* the camera P.
*/
inline bool isInFrontOfCamera(const Mat34 &P, const Vec4 &X){
inline bool isInFrontOfCamera(const Mat34 &P, const Vec4 &X) {
double condition_1 = P.row(2).dot(X) * X[3];
double condition_2 = X[2] * X[3];
if( condition_1 > 0 && condition_2 > 0 )
if (condition_1 > 0 && condition_2 > 0)
return true;
else
return false;
}
inline bool isInFrontOfCamera(const Mat34 &P, const Vec3 &X){
inline bool isInFrontOfCamera(const Mat34 &P, const Vec3 &X) {
Vec4 X_homo;
X_homo.segment<3>(0) = X;
X_homo(3) = 1;
@@ -200,14 +200,14 @@ inline bool isInFrontOfCamera(const Mat34 &P, const Vec3 &X){
* Transforms a 2D point from pixel image coordinates to a 2D point in
* normalized image coordinates.
*/
inline Vec2 ImageToNormImageCoordinates(Mat3 &Kinverse, Vec2 &x){
inline Vec2 ImageToNormImageCoordinates(Mat3 &Kinverse, Vec2 &x) {
Vec3 x_h = Kinverse*EuclideanToHomogeneous(x);
return HomogeneousToEuclidean( x_h );
}
/// Estimates the root mean square error (2D)
inline double RootMeanSquareError(const Mat2X &x_image,
const Mat4X &X_world,
inline double RootMeanSquareError(const Mat2X &x_image,
const Mat4X &X_world,
const Mat34 &P) {
size_t num_points = x_image.cols();
Mat2X dx = Project(P, X_world) - x_image;
@@ -215,10 +215,10 @@ inline double RootMeanSquareError(const Mat2X &x_image,
}
/// Estimates the root mean square error (2D)
inline double RootMeanSquareError(const Mat2X &x_image,
const Mat3X &X_world,
const Mat3 &K,
const Mat3 &R,
inline double RootMeanSquareError(const Mat2X &x_image,
const Mat3X &X_world,
const Mat3 &K,
const Mat3 &R,
const Vec3 &t) {
Mat34 P;
P_From_KRt(K, R, t, &P);
@@ -226,6 +226,6 @@ inline double RootMeanSquareError(const Mat2X &x_image,
Mat2X dx = Project(P, X_world) - x_image;
return dx.norm() / num_points;
}
} // namespace libmv
} // namespace libmv
#endif // LIBMV_MULTIVIEW_PROJECTION_H_

16
extern/libmv/libmv/multiview/resection.h vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2009 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -17,7 +17,7 @@
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
//
//
// Compute the projection matrix from a set of 3D points X and their
// projections x = PX in 2D. This is useful if a point cloud is reconstructed.
//
@@ -46,10 +46,10 @@ void Resection(const Matrix<T, 2, Dynamic> &x,
T xi = x(0, i);
T yi = x(1, i);
// See equation (7.2) on page 179 of H&Z.
design.template block<1,4>(2*i, 4) = -X.col(i).transpose();
design.template block<1,4>(2*i, 8) = yi*X.col(i).transpose();
design.template block<1,4>(2*i + 1, 0) = X.col(i).transpose();
design.template block<1,4>(2*i + 1, 8) = -xi*X.col(i).transpose();
design.template block<1, 4>(2*i, 4) = -X.col(i).transpose();
design.template block<1, 4>(2*i, 8) = yi*X.col(i).transpose();
design.template block<1, 4>(2*i + 1, 0) = X.col(i).transpose();
design.template block<1, 4>(2*i + 1, 8) = -xi*X.col(i).transpose();
}
Matrix<T, 12, 1> p;
Nullspace(&design, &p);

17
extern/libmv/libmv/multiview/triangulation.cc vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2007, 2008 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -18,9 +18,10 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/multiview/triangulation.h"
#include "libmv/numeric/numeric.h"
#include "libmv/multiview/projection.h"
#include "libmv/multiview/triangulation.h"
namespace libmv {
@@ -30,10 +31,10 @@ void TriangulateDLT(const Mat34 &P1, const Vec2 &x1,
Vec4 *X_homogeneous) {
Mat4 design;
for (int i = 0; i < 4; ++i) {
design(0,i) = x1(0) * P1(2,i) - P1(0,i);
design(1,i) = x1(1) * P1(2,i) - P1(1,i);
design(2,i) = x2(0) * P2(2,i) - P2(0,i);
design(3,i) = x2(1) * P2(2,i) - P2(1,i);
design(0, i) = x1(0) * P1(2, i) - P1(0, i);
design(1, i) = x1(1) * P1(2, i) - P1(1, i);
design(2, i) = x2(0) * P2(2, i) - P2(0, i);
design(3, i) = x2(1) * P2(2, i) - P2(1, i);
}
Nullspace(&design, X_homogeneous);
}

10
extern/libmv/libmv/multiview/triangulation.h vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2007, 2008 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -32,7 +32,7 @@ void TriangulateDLT(const Mat34 &P1, const Vec2 &x1,
void TriangulateDLT(const Mat34 &P1, const Vec2 &x1,
const Mat34 &P2, const Vec2 &x2,
Vec3 *X_euclidean);
} // namespace libmv
} // namespace libmv
#endif // LIBMV_MULTIVIEW_TRIANGULATION_H_

8
extern/libmv/libmv/numeric/dogleg.h vendored
View File

@@ -41,7 +41,7 @@ namespace libmv {
template<typename Function,
typename Jacobian = NumericJacobian<Function>,
typename Solver = Eigen::PartialPivLU<
Matrix<typename Function::FMatrixType::RealScalar,
Matrix<typename Function::FMatrixType::RealScalar,
Function::XMatrixType::RowsAtCompileTime,
Function::XMatrixType::RowsAtCompileTime> > >
class Dogleg {
@@ -52,8 +52,8 @@ class Dogleg {
typedef Matrix<typename Function::FMatrixType::RealScalar,
Function::FMatrixType::RowsAtCompileTime,
Function::XMatrixType::RowsAtCompileTime> JMatrixType;
typedef Matrix<typename JMatrixType::RealScalar,
JMatrixType::ColsAtCompileTime,
typedef Matrix<typename JMatrixType::RealScalar,
JMatrixType::ColsAtCompileTime,
JMatrixType::ColsAtCompileTime> AMatrixType;
enum Status {
@@ -75,7 +75,7 @@ class Dogleg {
: f_(f), df_(f) {}
struct SolverParameters {
SolverParameters()
SolverParameters()
: gradient_threshold(1e-16),
relative_step_threshold(1e-16),
error_threshold(1e-16),

4
extern/libmv/libmv/numeric/function_derivative.h vendored
View File

@@ -36,7 +36,7 @@ enum NumericJacobianMode {
FORWARD,
};
template<typename Function, NumericJacobianMode mode=CENTRAL>
template<typename Function, NumericJacobianMode mode = CENTRAL>
class NumericJacobian {
public:
typedef typename Function::XMatrixType Parameters;
@@ -58,7 +58,7 @@ class NumericJacobian {
XScalar mean_eps = eps.sum() / eps.rows();
if (mean_eps == XScalar(0)) {
// TODO(keir): Do something better here.
mean_eps = 1e-8; // ~sqrt(machine precision).
mean_eps = 1e-8; // ~sqrt(machine precision).
}
// TODO(keir): Elimininate this needless function evaluation for the
// central difference case.

14
extern/libmv/libmv/numeric/levenberg_marquardt.h vendored
View File

@@ -40,7 +40,7 @@ namespace libmv {
template<typename Function,
typename Jacobian = NumericJacobian<Function>,
typename Solver = Eigen::PartialPivLU<
Matrix<typename Function::FMatrixType::RealScalar,
Matrix<typename Function::FMatrixType::RealScalar,
Function::XMatrixType::RowsAtCompileTime,
Function::XMatrixType::RowsAtCompileTime> > >
class LevenbergMarquardt {
@@ -51,8 +51,8 @@ class LevenbergMarquardt {
typedef Matrix<typename Function::FMatrixType::RealScalar,
Function::FMatrixType::RowsAtCompileTime,
Function::XMatrixType::RowsAtCompileTime> JMatrixType;
typedef Matrix<typename JMatrixType::RealScalar,
JMatrixType::ColsAtCompileTime,
typedef Matrix<typename JMatrixType::RealScalar,
JMatrixType::ColsAtCompileTime,
JMatrixType::ColsAtCompileTime> AMatrixType;
// TODO(keir): Some of these knobs can be derived from each other and
@@ -69,7 +69,7 @@ class LevenbergMarquardt {
: f_(f), df_(f) {}
struct SolverParameters {
SolverParameters()
SolverParameters()
: gradient_threshold(1e-16),
relative_step_threshold(1e-16),
error_threshold(1e-16),
@@ -140,7 +140,7 @@ class LevenbergMarquardt {
if (solved && dx.norm() <= params.relative_step_threshold * x.norm()) {
results.status = RELATIVE_STEP_SIZE_TOO_SMALL;
break;
}
}
if (solved) {
x_new = x + dx;
// Rho is the ratio of the actual reduction in error to the reduction
@@ -156,8 +156,8 @@ class LevenbergMarquardt {
u = u*std::max(1/3., 1 - (tmp*tmp*tmp));
v = 2;
continue;
}
}
}
}
// Reject the update because either the normal equations failed to solve
// or the local linear model was not good (rho < 0). Instead, increase u
// to move closer to gradient descent.

6
extern/libmv/libmv/numeric/numeric.cc vendored
View File

@@ -121,15 +121,15 @@ void HorizontalStack(const Mat &left, const Mat &right, Mat *stacked) {
void MatrixColumn(const Mat &A, int i, Vec2 *v) {
assert(A.rows() == 2);
*v << A(0,i), A(1,i);
*v << A(0, i), A(1, i);
}
void MatrixColumn(const Mat &A, int i, Vec3 *v) {
assert(A.rows() == 3);
*v << A(0,i), A(1,i), A(2,i);
*v << A(0, i), A(1, i), A(2, i);
}
void MatrixColumn(const Mat &A, int i, Vec4 *v) {
assert(A.rows() == 4);
*v << A(0,i), A(1,i), A(2,i), A(3,i);
*v << A(0, i), A(1, i), A(2, i), A(3, i);
}
} // namespace libmv

69
extern/libmv/libmv/numeric/numeric.h vendored
View File

@@ -34,11 +34,11 @@
#include <Eigen/SVD>
#if (defined(_WIN32) || defined(__APPLE__) || defined(__FreeBSD__)) && !defined(__MINGW64__)
void static sincos (double x, double *sinx, double *cosx) {
static void sincos(double x, double *sinx, double *cosx) {
*sinx = sin(x);
*cosx = cos(x);
}
#endif //_WIN32 || __APPLE__
#endif // _WIN32 || __APPLE__
#if (defined(WIN32) || defined(WIN64)) && !defined(__MINGW32__)
inline long lround(double d) {
@@ -48,7 +48,7 @@
return (d>0) ? int(d+0.5) : int(d-0.5);
}
typedef unsigned int uint;
#endif //_WIN32
#endif // _WIN32
namespace libmv {
@@ -86,16 +86,16 @@ typedef Eigen::Matrix<double, 4, 4, Eigen::RowMajor> RMat4;
typedef Eigen::Matrix<double, 2, Eigen::Dynamic> Mat2X;
typedef Eigen::Matrix<double, 3, Eigen::Dynamic> Mat3X;
typedef Eigen::Matrix<double, 4, Eigen::Dynamic> Mat4X;
typedef Eigen::Matrix<double, Eigen::Dynamic, 2> MatX2;
typedef Eigen::Matrix<double, Eigen::Dynamic, 3> MatX3;
typedef Eigen::Matrix<double, Eigen::Dynamic, 4> MatX4;
typedef Eigen::Matrix<double, Eigen::Dynamic, 5> MatX5;
typedef Eigen::Matrix<double, Eigen::Dynamic, 6> MatX6;
typedef Eigen::Matrix<double, Eigen::Dynamic, 7> MatX7;
typedef Eigen::Matrix<double, Eigen::Dynamic, 8> MatX8;
typedef Eigen::Matrix<double, Eigen::Dynamic, 9> MatX9;
typedef Eigen::Matrix<double, Eigen::Dynamic,15> MatX15;
typedef Eigen::Matrix<double, Eigen::Dynamic,16> MatX16;
typedef Eigen::Matrix<double, Eigen::Dynamic, 2> MatX2;
typedef Eigen::Matrix<double, Eigen::Dynamic, 3> MatX3;
typedef Eigen::Matrix<double, Eigen::Dynamic, 4> MatX4;
typedef Eigen::Matrix<double, Eigen::Dynamic, 5> MatX5;
typedef Eigen::Matrix<double, Eigen::Dynamic, 6> MatX6;
typedef Eigen::Matrix<double, Eigen::Dynamic, 7> MatX7;
typedef Eigen::Matrix<double, Eigen::Dynamic, 8> MatX8;
typedef Eigen::Matrix<double, Eigen::Dynamic, 9> MatX9;
typedef Eigen::Matrix<double, Eigen::Dynamic, 15> MatX15;
typedef Eigen::Matrix<double, Eigen::Dynamic, 16> MatX16;
typedef Eigen::Vector2d Vec2;
typedef Eigen::Vector3d Vec3;
@@ -286,18 +286,21 @@ void MeanAndVarianceAlongRows(const Mat &A,
#if _WIN32
// TODO(bomboze): un-#if this for both platforms once tested under Windows
/* This solution was extensively discussed here http://forum.kde.org/viewtopic.php?f=74&t=61940 */
#define SUM_OR_DYNAMIC(x,y) (x==Eigen::Dynamic||y==Eigen::Dynamic)?Eigen::Dynamic:(x+y)
/* This solution was extensively discussed here
http://forum.kde.org/viewtopic.php?f=74&t=61940 */
#define SUM_OR_DYNAMIC(x, y) (x == Eigen::Dynamic || y == Eigen::Dynamic) ? Eigen::Dynamic : (x+y)
template<typename Derived1, typename Derived2>
struct hstack_return {
typedef typename Derived1::Scalar Scalar;
enum {
RowsAtCompileTime = Derived1::RowsAtCompileTime,
ColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::ColsAtCompileTime, Derived2::ColsAtCompileTime),
ColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::ColsAtCompileTime,
Derived2::ColsAtCompileTime),
Options = Derived1::Flags&Eigen::RowMajorBit ? Eigen::RowMajor : 0,
MaxRowsAtCompileTime = Derived1::MaxRowsAtCompileTime,
MaxColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxColsAtCompileTime, Derived2::MaxColsAtCompileTime)
MaxColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxColsAtCompileTime,
Derived2::MaxColsAtCompileTime)
};
typedef Eigen::Matrix<Scalar,
RowsAtCompileTime,
@@ -308,9 +311,10 @@ void MeanAndVarianceAlongRows(const Mat &A,
};
template<typename Derived1, typename Derived2>
typename hstack_return<Derived1,Derived2>::type
HStack (const Eigen::MatrixBase<Derived1>& lhs, const Eigen::MatrixBase<Derived2>& rhs) {
typename hstack_return<Derived1,Derived2>::type res;
typename hstack_return<Derived1, Derived2>::type
HStack(const Eigen::MatrixBase<Derived1>& lhs,
const Eigen::MatrixBase<Derived2>& rhs) {
typename hstack_return<Derived1, Derived2>::type res;
res.resize(lhs.rows(), lhs.cols()+rhs.cols());
res << lhs, rhs;
return res;
@@ -321,10 +325,12 @@ void MeanAndVarianceAlongRows(const Mat &A,
struct vstack_return {
typedef typename Derived1::Scalar Scalar;
enum {
RowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::RowsAtCompileTime, Derived2::RowsAtCompileTime),
RowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::RowsAtCompileTime,
Derived2::RowsAtCompileTime),
ColsAtCompileTime = Derived1::ColsAtCompileTime,
Options = Derived1::Flags&Eigen::RowMajorBit ? Eigen::RowMajor : 0,
MaxRowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxRowsAtCompileTime, Derived2::MaxRowsAtCompileTime),
MaxRowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxRowsAtCompileTime,
Derived2::MaxRowsAtCompileTime),
MaxColsAtCompileTime = Derived1::MaxColsAtCompileTime
};
typedef Eigen::Matrix<Scalar,
@@ -336,16 +342,17 @@ void MeanAndVarianceAlongRows(const Mat &A,
};
template<typename Derived1, typename Derived2>
typename vstack_return<Derived1,Derived2>::type
VStack (const Eigen::MatrixBase<Derived1>& lhs, const Eigen::MatrixBase<Derived2>& rhs) {
typename vstack_return<Derived1,Derived2>::type res;
typename vstack_return<Derived1, Derived2>::type
VStack(const Eigen::MatrixBase<Derived1>& lhs,
const Eigen::MatrixBase<Derived2>& rhs) {
typename vstack_return<Derived1, Derived2>::type res;
res.resize(lhs.rows()+rhs.rows(), lhs.cols());
res << lhs, rhs;
return res;
};
#else //_WIN32
#else // _WIN32
// Since it is not possible to typedef privately here, use a macro.
// Always take dynamic columns if either side is dynamic.
@@ -400,7 +407,7 @@ void MeanAndVarianceAlongRows(const Mat &A,
}
#undef COLS
#undef ROWS
#endif //_WIN32
#endif // _WIN32
@@ -454,7 +461,7 @@ inline bool isnan(double i) {
/// and negative values
template <typename FloatType>
FloatType ceil0(const FloatType& value) {
FloatType result = std::ceil( std::fabs( value ) );
FloatType result = std::ceil(std::fabs(value));
return (value < 0.0) ? -result : result;
}
@@ -470,8 +477,8 @@ inline Mat3 SkewMat(const Vec3 &x) {
/// the first two (independent) lines
inline Mat23 SkewMatMinimal(const Vec2 &x) {
Mat23 skew;
skew << 0,-1, x(1),
1, 0, -x(0);
skew << 0, -1, x(1),
1, 0, -x(0);
return skew;
}
@@ -485,6 +492,6 @@ inline Mat3 RotationFromEulerVector(Vec3 euler_vector) {
Mat3 w_hat = CrossProductMatrix(w);
return Mat3::Identity() + w_hat*sin(theta) + w_hat*w_hat*(1 - cos(theta));
}
} // namespace libmv
} // namespace libmv
#endif // LIBMV_NUMERIC_NUMERIC_H

6
extern/libmv/libmv/numeric/poly.cc vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2007, 2008 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

28
extern/libmv/libmv/numeric/poly.h vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2007, 2008 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -36,7 +36,7 @@ namespace libmv {
//
// The GSL cubic solver was used as a reference for this routine.
template<typename Real>
int SolveCubicPolynomial(Real a, Real b, Real c,
int SolveCubicPolynomial(Real a, Real b, Real c,
Real *x0, Real *x1, Real *x2) {
Real q = a * a - 3 * b;
Real r = 2 * a * a * a - 9 * a * b + 27 * c;
@@ -52,7 +52,7 @@ int SolveCubicPolynomial(Real a, Real b, Real c,
if (R == 0 && Q == 0) {
// Tripple root in one place.
*x0 = *x1 = *x2 = -a / 3 ;
*x0 = *x1 = *x2 = -a / 3;
return 3;
} else if (CR2 == CQ3) {
@@ -62,7 +62,7 @@ int SolveCubicPolynomial(Real a, Real b, Real c,
// Due to finite precision some double roots may be missed, and considered
// to be a pair of complex roots z = x +/- epsilon i close to the real
// axis.
Real sqrtQ = sqrt (Q);
Real sqrtQ = sqrt(Q);
if (R > 0) {
*x0 = -2 * sqrtQ - a / 3;
*x1 = sqrtQ - a / 3;
@@ -76,13 +76,13 @@ int SolveCubicPolynomial(Real a, Real b, Real c,
} else if (CR2 < CQ3) {
// This case is equivalent to R2 < Q3.
Real sqrtQ = sqrt (Q);
Real sqrtQ = sqrt(Q);
Real sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
Real theta = acos (R / sqrtQ3);
Real theta = acos(R / sqrtQ3);
Real norm = -2 * sqrtQ;
*x0 = norm * cos (theta / 3) - a / 3;
*x1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3;
*x2 = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3;
*x0 = norm * cos(theta / 3) - a / 3;
*x1 = norm * cos((theta + 2.0 * M_PI) / 3) - a / 3;
*x2 = norm * cos((theta - 2.0 * M_PI) / 3) - a / 3;
// Put the roots in ascending order.
if (*x0 > *x1) {
@@ -95,10 +95,10 @@ int SolveCubicPolynomial(Real a, Real b, Real c,
}
}
return 3;
}
}
Real sgnR = (R >= 0 ? 1 : -1);
Real A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0/3.0);
Real B = Q / A ;
Real A = -sgnR * pow(fabs(R) + sqrt(R2 - Q3), 1.0/3.0);
Real B = Q / A;
*x0 = A + B - a / 3;
return 1;
}

15
extern/libmv/libmv/simple_pipeline/bundle.cc vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2011, 2012, 2013 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -18,6 +18,8 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/simple_pipeline/bundle.h"
#include <map>
#include "ceres/ceres.h"
@@ -29,7 +31,6 @@
#include "libmv/multiview/projection.h"
#include "libmv/numeric/numeric.h"
#include "libmv/simple_pipeline/camera_intrinsics.h"
#include "libmv/simple_pipeline/bundle.h"
#include "libmv/simple_pipeline/reconstruction.h"
#include "libmv/simple_pipeline/tracks.h"
@@ -125,8 +126,8 @@ struct OpenCVReprojectionError {
// single parameter block.
struct RotationMatrixPlus {
template<typename T>
bool operator()(const T* R_array, // Rotation 3x3 col-major.
const T* delta, // Angle-axis delta
bool operator()(const T* R_array, // Rotation 3x3 col-major.
const T* delta, // Angle-axis delta
T* R_plus_delta_array) const {
T angle_axis[3];
@@ -293,7 +294,7 @@ void EuclideanBundleCommonIntrinsics(const Tracks &tracks,
}
LG << "Number of residuals: " << num_residuals;
if(!num_residuals) {
if (!num_residuals) {
LG << "Skipping running minimizer with zero residuals";
return;
}

13
extern/libmv/libmv/simple_pipeline/bundle.h vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2011 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -89,14 +89,15 @@ enum BundleIntrinsics {
BUNDLE_TANGENTIAL = 48,
};
enum BundleConstraints {
BUNDLE_NO_CONSTRAINTS = 0,
BUNDLE_NO_TRANSLATION = 1,
BUNDLE_NO_CONSTRAINTS = 0,
BUNDLE_NO_TRANSLATION = 1,
};
void EuclideanBundleCommonIntrinsics(const Tracks &tracks,
int bundle_intrinsics,
EuclideanReconstruction *reconstruction,
CameraIntrinsics *intrinsics,
int bundle_constraints = BUNDLE_NO_CONSTRAINTS);
int bundle_constraints =
BUNDLE_NO_CONSTRAINTS);
/*!
Refine camera poses and 3D coordinates using bundle adjustment.

29
extern/libmv/libmv/simple_pipeline/callbacks.cc vendored
View File

@@ -1,29 +0,0 @@
// Copyright (c) 2011 libmv authors.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/simple_pipeline/callbacks.h"
namespace libmv {
void ProgressUpdateCallback::invoke(double progress, const char* message)
{
}
} // namespace libmv

3
extern/libmv/libmv/simple_pipeline/callbacks.h vendored
View File

@@ -25,7 +25,8 @@ namespace libmv {
class ProgressUpdateCallback {
public:
virtual void invoke(double progress, const char *message);
virtual ~ProgressUpdateCallback() {}
virtual void invoke(double progress, const char *message) = 0;
};
} // namespace libmv

151
extern/libmv/libmv/simple_pipeline/camera_intrinsics.cc vendored
View File

@@ -25,7 +25,7 @@ namespace libmv {
struct Offset {
short ix, iy;
unsigned char fx,fy;
unsigned char fx, fy;
};
struct Grid {
@@ -34,8 +34,7 @@ struct Grid {
double overscan;
};
static struct Grid *copyGrid(struct Grid *from)
{
static struct Grid *copyGrid(struct Grid *from) {
struct Grid *to = NULL;
if (from) {
@@ -74,8 +73,7 @@ CameraIntrinsics::CameraIntrinsics(const CameraIntrinsics &from)
k3_(from.k3_),
p1_(from.p1_),
p2_(from.p2_),
threads_(from.threads_)
{
threads_(from.threads_) {
distort_ = copyGrid(from.distort_);
undistort_ = copyGrid(from.undistort_);
}
@@ -122,9 +120,8 @@ void CameraIntrinsics::SetTangentialDistortion(double p1, double p2) {
FreeLookupGrid();
}
void CameraIntrinsics::SetThreads(int threads)
{
threads_ = threads;
void CameraIntrinsics::SetThreads(int threads) {
threads_ = threads;
}
void CameraIntrinsics::ApplyIntrinsics(double normalized_x,
@@ -189,68 +186,76 @@ void CameraIntrinsics::InvertIntrinsics(double image_x,
// TODO(MatthiasF): downsample lookup
template<typename WarpFunction>
void CameraIntrinsics::ComputeLookupGrid(Grid* grid, int width, int height, double overscan) {
void CameraIntrinsics::ComputeLookupGrid(Grid* grid, int width, int height,
double overscan) {
double w = (double)width / (1 + overscan);
double h = (double)height / (1 + overscan);
double aspx = (double)w / image_width_;
double aspy = (double)h / image_height_;
#if defined(_OPENMP)
#pragma omp parallel for schedule(dynamic) num_threads(threads_) if (threads_ > 1 && height > 100)
#endif
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
double src_x = (x - 0.5 * overscan * w) / aspx, src_y = (y - 0.5 * overscan * h) / aspy;
double src_x = (x - 0.5 * overscan * w) / aspx,
src_y = (y - 0.5 * overscan * h) / aspy;
double warp_x, warp_y;
WarpFunction(this,src_x,src_y,&warp_x,&warp_y);
WarpFunction(this, src_x, src_y, &warp_x, &warp_y);
warp_x = warp_x*aspx + 0.5 * overscan * w;
warp_y = warp_y*aspy + 0.5 * overscan * h;
int ix = int(warp_x), iy = int(warp_y);
int fx = round((warp_x-ix)*256), fy = round((warp_y-iy)*256);
if(fx == 256) { fx=0; ix++; }
if(fy == 256) { fy=0; iy++; }
if (fx == 256) { fx = 0; ix++; } // NOLINT
if (fy == 256) { fy = 0; iy++; } // NOLINT
// Use nearest border pixel
if( ix < 0 ) { ix = 0, fx = 0; }
if( iy < 0 ) { iy = 0, fy = 0; }
if( ix >= width-2 ) ix = width-2;
if( iy >= height-2 ) iy = height-2;
if (ix < 0) { ix = 0, fx = 0; } // NOLINT
if (iy < 0) { iy = 0, fy = 0; } // NOLINT
if (ix >= width - 2) ix = width-2;
if (iy >= height - 2) iy = height-2;
Offset offset = { (short)(ix-x), (short)(iy-y), (unsigned char)fx, (unsigned char)fy };
Offset offset = { (short)(ix-x), (short)(iy-y),
(unsigned char)fx, (unsigned char)fy };
grid->offset[y*width+x] = offset;
}
}
}
// TODO(MatthiasF): cubic B-Spline image sampling, bilinear lookup
template<typename T,int N>
template<typename T, int N>
static void Warp(const Grid* grid, const T* src, T* dst,
int width, int height, int threads) {
(void) threads; // Ignored if OpenMP is disabled
#if defined(_OPENMP)
#pragma omp parallel for schedule(dynamic) num_threads(threads) if (threads > 1 && height > 100)
#endif
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
Offset offset = grid->offset[y*width+x];
const T* s = &src[((y+offset.iy)*width+(x+offset.ix))*N];
for (int i = 0; i < N; i++) {
dst[(y*width+x)*N+i] = ((s[ i] * (256-offset.fx) + s[ N+i] * offset.fx) * (256-offset.fy)
+(s[width*N+i] * (256-offset.fx) + s[width*N+N+i] * offset.fx) * offset.fy) / (256*256);
dst[(y*width+x)*N+i] = ((s[ i] * (256-offset.fx) + s[ N+i] * offset.fx) * (256-offset.fy) // NOLINT
+(s[width*N+i] * (256-offset.fx) + s[width*N+N+i] * offset.fx) * offset.fy) / (256*256); // NOLINT
}
}
}
}
void CameraIntrinsics::FreeLookupGrid() {
if(distort_) {
if (distort_) {
delete distort_->offset;
delete distort_;
distort_ = NULL;
}
if(undistort_) {
if (undistort_) {
delete undistort_->offset;
delete undistort_;
undistort_ = NULL;
}
}
// FIXME: C++ templates limitations makes thing complicated, but maybe there is a simpler method.
// FIXME: C++ templates limitations makes thing complicated,
// but maybe there is a simpler method.
struct ApplyIntrinsicsFunction {
ApplyIntrinsicsFunction(CameraIntrinsics* intrinsics, double x, double y,
double *warp_x, double *warp_y) {
@@ -263,16 +268,21 @@ struct ApplyIntrinsicsFunction {
struct InvertIntrinsicsFunction {
InvertIntrinsicsFunction(CameraIntrinsics* intrinsics, double x, double y,
double *warp_x, double *warp_y) {
intrinsics->InvertIntrinsics(x,y,warp_x,warp_y);
*warp_x = *warp_x*intrinsics->focal_length_x()+intrinsics->principal_point_x();
*warp_y = *warp_y*intrinsics->focal_length_y()+intrinsics->principal_point_y();
intrinsics->InvertIntrinsics(x, y, warp_x, warp_y);
*warp_x = *warp_x * intrinsics->focal_length_x() +
intrinsics->principal_point_x();
*warp_y = *warp_y * intrinsics->focal_length_y() +
intrinsics->principal_point_y();
}
};
void CameraIntrinsics::CheckDistortLookupGrid(int width, int height, double overscan)
{
if(distort_) {
if(distort_->width != width || distort_->height != height || distort_->overscan != overscan) {
void CameraIntrinsics::CheckDistortLookupGrid(int width, int height,
double overscan) {
if (distort_) {
if (distort_->width != width || distort_->height != height ||
distort_->overscan != overscan) {
delete [] distort_->offset;
distort_->offset = NULL;
}
@@ -281,9 +291,10 @@ void CameraIntrinsics::CheckDistortLookupGrid(int width, int height, double over
distort_->offset = NULL;
}
if(!distort_->offset) {
distort_->offset = new Offset[width*height];
ComputeLookupGrid<InvertIntrinsicsFunction>(distort_,width,height,overscan);
if (!distort_->offset) {
distort_->offset = new Offset[width * height];
ComputeLookupGrid<InvertIntrinsicsFunction>(distort_, width,
height, overscan);
}
distort_->width = width;
@@ -291,10 +302,11 @@ void CameraIntrinsics::CheckDistortLookupGrid(int width, int height, double over
distort_->overscan = overscan;
}
void CameraIntrinsics::CheckUndistortLookupGrid(int width, int height, double overscan)
{
if(undistort_) {
if(undistort_->width != width || undistort_->height != height || undistort_->overscan != overscan) {
void CameraIntrinsics::CheckUndistortLookupGrid(int width, int height,
double overscan) {
if (undistort_) {
if (undistort_->width != width || undistort_->height != height ||
undistort_->overscan != overscan) {
delete [] undistort_->offset;
undistort_->offset = NULL;
}
@@ -303,9 +315,10 @@ void CameraIntrinsics::CheckUndistortLookupGrid(int width, int height, double ov
undistort_->offset = NULL;
}
if(!undistort_->offset) {
undistort_->offset = new Offset[width*height];
ComputeLookupGrid<ApplyIntrinsicsFunction>(undistort_,width,height,overscan);
if (!undistort_->offset) {
undistort_->offset = new Offset[width * height];
ComputeLookupGrid<ApplyIntrinsicsFunction>(undistort_, width,
height, overscan);
}
undistort_->width = width;
@@ -313,39 +326,53 @@ void CameraIntrinsics::CheckUndistortLookupGrid(int width, int height, double ov
undistort_->overscan = overscan;
}
void CameraIntrinsics::Distort(const float* src, float* dst, int width, int height, double overscan, int channels) {
void CameraIntrinsics::Distort(const float* src, float* dst,
int width, int height,
double overscan,
int channels) {
CheckDistortLookupGrid(width, height, overscan);
if(channels==1) Warp<float,1>(distort_,src,dst,width,height,threads_);
else if(channels==2) Warp<float,2>(distort_,src,dst,width,height,threads_);
else if(channels==3) Warp<float,3>(distort_,src,dst,width,height,threads_);
else if(channels==4) Warp<float,4>(distort_,src,dst,width,height,threads_);
if (channels==1) Warp<float,1>(distort_, src, dst, width, height, threads_); // NOLINT
else if (channels==2) Warp<float,2>(distort_, src, dst, width, height, threads_); // NOLINT
else if (channels==3) Warp<float,3>(distort_, src, dst, width, height, threads_); // NOLINT
else if (channels==4) Warp<float,4>(distort_, src, dst, width, height, threads_); // NOLINT
//else assert("channels must be between 1 and 4");
}
void CameraIntrinsics::Distort(const unsigned char* src, unsigned char* dst, int width, int height, double overscan, int channels) {
void CameraIntrinsics::Distort(const unsigned char* src,
unsigned char* dst,
int width, int height,
double overscan,
int channels) {
CheckDistortLookupGrid(width, height, overscan);
if(channels==1) Warp<unsigned char,1>(distort_,src,dst,width,height,threads_);
else if(channels==2) Warp<unsigned char,2>(distort_,src,dst,width,height,threads_);
else if(channels==3) Warp<unsigned char,3>(distort_,src,dst,width,height,threads_);
else if(channels==4) Warp<unsigned char,4>(distort_,src,dst,width,height,threads_);
if (channels == 1) Warp<unsigned char,1>(distort_, src, dst, width, height, threads_); // NOLINT
else if (channels == 2) Warp<unsigned char,2>(distort_, src, dst, width, height, threads_); // NOLINT
else if (channels == 3) Warp<unsigned char,3>(distort_, src, dst, width, height, threads_); // NOLINT
else if (channels == 4) Warp<unsigned char,4>(distort_, src, dst, width, height, threads_); // NOLINT
//else assert("channels must be between 1 and 4");
}
void CameraIntrinsics::Undistort(const float* src, float* dst, int width, int height, double overscan, int channels) {
void CameraIntrinsics::Undistort(const float* src, float* dst,
int width, int height,
double overscan,
int channels) {
CheckUndistortLookupGrid(width, height, overscan);
if(channels==1) Warp<float,1>(undistort_,src,dst,width,height,threads_);
else if(channels==2) Warp<float,2>(undistort_,src,dst,width,height,threads_);
else if(channels==3) Warp<float,3>(undistort_,src,dst,width,height,threads_);
else if(channels==4) Warp<float,4>(undistort_,src,dst,width,height,threads_);
if (channels == 1) Warp<float,1>(undistort_, src, dst, width, height, threads_); // NOLINT
else if (channels == 2) Warp<float,2>(undistort_, src, dst, width, height, threads_); // NOLINT
else if (channels == 3) Warp<float,3>(undistort_, src, dst, width, height, threads_); // NOLINT
else if (channels == 4) Warp<float,4>(undistort_, src, dst, width, height, threads_); // NOLINT
//else assert("channels must be between 1 and 4");
}
void CameraIntrinsics::Undistort(const unsigned char* src, unsigned char* dst, int width, int height, double overscan, int channels) {
void CameraIntrinsics::Undistort(const unsigned char* src,
unsigned char* dst,
int width, int height,
double overscan,
int channels) {
CheckUndistortLookupGrid(width, height, overscan);
if(channels==1) Warp<unsigned char,1>(undistort_,src,dst,width,height,threads_);
else if(channels==2) Warp<unsigned char,2>(undistort_,src,dst,width,height,threads_);
else if(channels==3) Warp<unsigned char,3>(undistort_,src,dst,width,height,threads_);
else if(channels==4) Warp<unsigned char,4>(undistort_,src,dst,width,height,threads_);
if (channels == 1) Warp<unsigned char,1>(undistort_, src, dst, width, height, threads_); // NOLINT
else if (channels == 2) Warp<unsigned char,2>(undistort_, src, dst, width, height, threads_); // NOLINT
else if (channels == 3) Warp<unsigned char,3>(undistort_, src, dst, width, height, threads_); // NOLINT
else if (channels == 4) Warp<unsigned char,4>(undistort_, src, dst, width, height, threads_); // NOLINT
//else assert("channels must be between 1 and 4");
}

5
extern/libmv/libmv/simple_pipeline/camera_intrinsics.h vendored
View File

@@ -135,7 +135,10 @@ class CameraIntrinsics {
int width, int height, double overscan, int channels);
private:
template<typename WarpFunction> void ComputeLookupGrid(struct Grid* grid, int width, int height, double overscan);
template<typename WarpFunction> void ComputeLookupGrid(struct Grid* grid,
int width,
int height,
double overscan);
void CheckUndistortLookupGrid(int width, int height, double overscan);
void CheckDistortLookupGrid(int width, int height, double overscan);
void FreeLookupGrid();

121
extern/libmv/libmv/simple_pipeline/detect.cc vendored
View File

@@ -35,37 +35,40 @@ namespace libmv {
typedef unsigned int uint;
static int featurecmp(const void *a_v, const void *b_v)
{
static int featurecmp(const void *a_v, const void *b_v) {
Feature *a = (Feature*)a_v;
Feature *b = (Feature*)b_v;
return b->score - a->score;
}
std::vector<Feature> DetectFAST(const unsigned char* data, int width, int height, int stride,
int min_trackness, int min_distance) {
std::vector<Feature> DetectFAST(const unsigned char* data,
int width, int height,
int stride,
int min_trackness,
int min_distance) {
std::vector<Feature> features;
// TODO(MatthiasF): Support targetting a feature count (binary search trackness)
int num_features;
xy* all = fast9_detect(data, width, height,
stride, min_trackness, &num_features);
if(num_features == 0) {
if (num_features == 0) {
free(all);
return features;
}
int* scores = fast9_score(data, stride, all, num_features, min_trackness);
// TODO: merge with close feature suppression
// TODO(MatthiasF): merge with close feature suppression
xy* nonmax = nonmax_suppression(all, scores, num_features, &num_features);
free(all);
// Remove too close features
// TODO(MatthiasF): A resolution independent parameter would be better than distance
// e.g. a coefficient going from 0 (no minimal distance) to 1 (optimal circle packing)
// TODO(MatthiasF): A resolution independent parameter would be better than
// distance e.g. a coefficient going from 0 (no minimal distance) to 1
// (optimal circle packing)
// FIXME(MatthiasF): this method will not necessarily give all maximum markers
if(num_features) {
if (num_features) {
Feature *all_features = new Feature[num_features];
for(int i = 0; i < num_features; ++i) {
for (int i = 0; i < num_features; ++i) {
Feature a = { (float)nonmax[i].x, (float)nonmax[i].y, (float)scores[i], 0 };
all_features[i] = a;
}
@@ -75,15 +78,15 @@ std::vector<Feature> DetectFAST(const unsigned char* data, int width, int height
features.reserve(num_features);
int prev_score = all_features[0].score;
for(int i = 0; i < num_features; ++i) {
for (int i = 0; i < num_features; ++i) {
bool ok = true;
Feature a = all_features[i];
if(a.score>prev_score)
if (a.score>prev_score)
abort();
prev_score = a.score;
// compare each feature against filtered set
for(int j = 0; j < features.size(); j++) {
for (int j = 0; j < features.size(); j++) {
Feature& b = features[j];
if ( (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y) < min_distance*min_distance ) {
// already a nearby feature
@@ -92,7 +95,7 @@ std::vector<Feature> DetectFAST(const unsigned char* data, int width, int height
}
}
if(ok) {
if (ok) {
// add the new feature
features.push_back(a);
}
@@ -106,19 +109,22 @@ std::vector<Feature> DetectFAST(const unsigned char* data, int width, int height
}
#ifdef __SSE2__
static uint SAD(const ubyte* imageA, const ubyte* imageB, int strideA, int strideB) {
static uint SAD(const ubyte* imageA, const ubyte* imageB,
int strideA, int strideB) {
__m128i a = _mm_setzero_si128();
for(int i = 0; i < 16; i++) {
a = _mm_adds_epu16(a, _mm_sad_epu8( _mm_loadu_si128((__m128i*)(imageA+i*strideA)),
_mm_loadu_si128((__m128i*)(imageB+i*strideB))));
for (int i = 0; i < 16; i++) {
a = _mm_adds_epu16(a,
_mm_sad_epu8(_mm_loadu_si128((__m128i*)(imageA+i*strideA)),
_mm_loadu_si128((__m128i*)(imageB+i*strideB))));
}
return _mm_extract_epi16(a,0) + _mm_extract_epi16(a,4);
return _mm_extract_epi16(a, 0) + _mm_extract_epi16(a, 4);
}
#else
static uint SAD(const ubyte* imageA, const ubyte* imageB, int strideA, int strideB) {
uint sad=0;
for(int i = 0; i < 16; i++) {
for(int j = 0; j < 16; j++) {
static uint SAD(const ubyte* imageA, const ubyte* imageB,
int strideA, int strideB) {
uint sad = 0;
for (int i = 0; i < 16; i++) {
for (int j = 0; j < 16; j++) {
sad += abs((int)imageA[i*strideA+j] - imageB[i*strideB+j]);
}
}
@@ -126,59 +132,68 @@ static uint SAD(const ubyte* imageA, const ubyte* imageB, int strideA, int strid
}
#endif
void DetectMORAVEC(ubyte* image, int stride, int width, int height, Feature* detected, int* count, int distance, ubyte* pattern) {
void DetectMORAVEC(ubyte* image,
int stride, int width, int height,
Feature* detected, int* count,
int distance,
ubyte* pattern) {
unsigned short histogram[256];
memset(histogram,0,sizeof(histogram));
memset(histogram, 0, sizeof(histogram));
ubyte* scores = new ubyte[width*height];
memset(scores,0,width*height);
const int r = 1; //radius for self similarity comparison
for(int y=distance; y<height-distance; y++) {
for(int x=distance; x<width-distance; x++) {
memset(scores, 0, width*height);
const int r = 1; // radius for self similarity comparison
for (int y = distance; y < height-distance; y++) {
for (int x = distance; x < width-distance; x++) {
ubyte* s = &image[y*stride+x];
int score = // low self-similarity with overlapping patterns //OPTI: load pattern once
int score = // low self-similarity with overlapping patterns
// OPTI: load pattern once
SAD(s, s-r*stride-r, stride, stride)+SAD(s, s-r*stride, stride, stride)+SAD(s, s-r*stride+r, stride, stride)+
SAD(s, s -r, stride, stride)+ SAD(s, s +r, stride, stride)+
SAD(s, s+r*stride-r, stride, stride)+SAD(s, s+r*stride, stride, stride)+SAD(s, s+r*stride+r, stride, stride);
score /= 256; // normalize
if(pattern) score -= SAD(s, pattern, stride, 16); // find only features similar to pattern
if(score<=16) continue; // filter very self-similar features
score -= 16; // translate to score/histogram values
if(score>255) score=255; // clip
score /= 256; // normalize
if (pattern) // find only features similar to pattern
score -= SAD(s, pattern, stride, 16);
if (score <= 16) continue; // filter very self-similar features
score -= 16; // translate to score/histogram values
if (score>255) score=255; // clip
ubyte* c = &scores[y*width+x];
for(int i=-distance; i<0; i++) {
for(int j=-distance; j<distance; j++) {
for (int i = -distance; i < 0; i++) {
for (int j = -distance; j < distance; j++) {
int s = c[i*width+j];
if(s == 0) continue;
if(s >= score) goto nonmax;
c[i*width+j]=0, histogram[s]--;
if (s == 0) continue;
if (s >= score) goto nonmax;
c[i*width+j] = 0;
histogram[s]--;
}
}
for(int i=0, j=-distance; j<0; j++) {
for (int i = 0, j = -distance; j < 0; j++) {
int s = c[i*width+j];
if(s == 0) continue;
if(s >= score) goto nonmax;
c[i*width+j]=0, histogram[s]--;
if (s == 0) continue;
if (s >= score) goto nonmax;
c[i*width+j] = 0;
histogram[s]--;
}
c[0] = score, histogram[score]++;
nonmax:;
nonmax:
{ } // Do nothing.
}
}
int min=255, total=0;
for(; min>0; min--) {
int min = 255, total = 0;
for (; min > 0; min--) {
int h = histogram[min];
if(total+h > *count) break;
if (total+h > *count) break;
total += h;
}
int i=0;
for(int y=16; y<height-16; y++) {
for(int x=16; x<width-16; x++) {
int i = 0;
for (int y = 16; y < height-16; y++) {
for (int x = 16; x < width-16; x++) {
int s = scores[y*width+x];
Feature f = { (float)x+8.0f, (float)y+8.0f, (float)s, 16 };
if(s>min) detected[i++] = f;
if (s > min) detected[i++] = f;
}
}
*count = i;
delete[] scores;
}
}
} // namespace libmv

6
extern/libmv/libmv/simple_pipeline/detect.h vendored
View File

@@ -88,8 +88,10 @@ std::vector<Feature> DetectFAST(const unsigned char* data, int width, int height
\note \a You can crop the image (to avoid detecting markers near the borders) without copying:
image += marginY*stride+marginX, width -= 2*marginX, height -= 2*marginY;
*/
void DetectMORAVEC(ubyte* image, int stride, int width, int height, Feature* detected, int* count, int distance /*=32*/, ubyte* pattern /*=0*/);
void DetectMORAVEC(ubyte* image, int stride, int width, int height,
Feature* detected, int* count, int distance /*=32*/,
ubyte* pattern /*=0*/);
}
} // namespace libmv
#endif

14
extern/libmv/libmv/simple_pipeline/initialize_reconstruction.cc vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2011 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -18,13 +18,14 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/simple_pipeline/initialize_reconstruction.h"
#include "libmv/base/vector.h"
#include "libmv/logging/logging.h"
#include "libmv/multiview/fundamental.h"
#include "libmv/multiview/projection.h"
#include "libmv/numeric/levenberg_marquardt.h"
#include "libmv/numeric/numeric.h"
#include "libmv/simple_pipeline/initialize_reconstruction.h"
#include "libmv/simple_pipeline/reconstruction.h"
#include "libmv/simple_pipeline/tracks.h"
@@ -109,7 +110,8 @@ namespace {
Mat3 DecodeF(const Vec9 &encoded_F) {
// Decode F and force it to be rank 2.
Map<const Mat3> full_rank_F(encoded_F.data(), 3, 3);
Eigen::JacobiSVD<Mat3> svd(full_rank_F, Eigen::ComputeFullU | Eigen::ComputeFullV);
Eigen::JacobiSVD<Mat3> svd(full_rank_F,
Eigen::ComputeFullU | Eigen::ComputeFullV);
Vec3 diagonal = svd.singularValues();
diagonal(2) = 0;
Mat3 F = svd.matrixU() * diagonal.asDiagonal() * svd.matrixV().transpose();
@@ -130,7 +132,7 @@ struct FundamentalSampsonCostFunction {
Vec operator()(const Vec9 &encoded_F) const {
// Decode F and force it to be rank 2.
Mat3 F = DecodeF(encoded_F);
Vec residuals(markers.size() / 2);
residuals.setZero();
for (int i = 0; i < markers.size() / 2; ++i) {

6
extern/libmv/libmv/simple_pipeline/initialize_reconstruction.h vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2011 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

7
extern/libmv/libmv/simple_pipeline/intersect.cc vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2011 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -25,7 +25,6 @@
#include "libmv/multiview/projection.h"
#include "libmv/multiview/triangulation.h"
#include "libmv/multiview/nviewtriangulation.h"
#include "libmv/multiview/projection.h"
#include "libmv/numeric/numeric.h"
#include "libmv/numeric/levenberg_marquardt.h"
#include "libmv/simple_pipeline/reconstruction.h"

6
extern/libmv/libmv/simple_pipeline/intersect.h vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2011 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

12
extern/libmv/libmv/simple_pipeline/modal_solver.cc vendored
View File

@@ -34,7 +34,7 @@
namespace libmv {
namespace {
void ProjectMarkerOnSphere(Marker &marker, Vec3 &X) {
void ProjectMarkerOnSphere(const Marker &marker, Vec3 &X) {
X(0) = marker.x;
X(1) = marker.y;
X(2) = 1.0;
@@ -43,8 +43,7 @@ void ProjectMarkerOnSphere(Marker &marker, Vec3 &X) {
}
void ModalSolverLogProress(ProgressUpdateCallback *update_callback,
double progress)
{
double progress) {
if (update_callback) {
char message[256];
@@ -56,13 +55,14 @@ void ModalSolverLogProress(ProgressUpdateCallback *update_callback,
}
struct ModalReprojectionError {
ModalReprojectionError(double observed_x, double observed_y, Vec3 &bundle)
ModalReprojectionError(double observed_x,
double observed_y,
const Vec3 &bundle)
: observed_x(observed_x), observed_y(observed_y), bundle(bundle) { }
template <typename T>
bool operator()(const T* quaternion, // Rotation quaternion
T* residuals) const {
T R[9];
ceres::QuaternionToRotation(quaternion, R);
@@ -96,7 +96,7 @@ struct ModalReprojectionError {
};
} // namespace
void ModalSolver(Tracks &tracks,
void ModalSolver(const Tracks &tracks,
EuclideanReconstruction *reconstruction,
ProgressUpdateCallback *update_callback) {
int max_image = tracks.MaxImage();

2
extern/libmv/libmv/simple_pipeline/modal_solver.h vendored
View File

@@ -39,7 +39,7 @@ namespace libmv {
Reconstructed cameras and projected bundles would be added to reconstruction
object.
*/
void ModalSolver(Tracks &tracks,
void ModalSolver(const Tracks &tracks,
EuclideanReconstruction *reconstruction,
ProgressUpdateCallback *update_callback = NULL);

35
extern/libmv/libmv/simple_pipeline/pipeline.cc vendored
View File

@@ -18,10 +18,11 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/simple_pipeline/pipeline.h"
#include <cstdio>
#include "libmv/logging/logging.h"
#include "libmv/simple_pipeline/pipeline.h"
#include "libmv/simple_pipeline/bundle.h"
#include "libmv/simple_pipeline/intersect.h"
#include "libmv/simple_pipeline/resect.h"
@@ -92,6 +93,9 @@ struct ProjectivePipelineRoutines {
static bool Resect(const ReconstructionOptions &options,
const vector<Marker> &markers,
ProjectiveReconstruction *reconstruction, bool final_pass) {
(void) options; // Ignored.
(void) final_pass; // Ignored.
return ProjectiveResect(markers, reconstruction);
}
@@ -120,17 +124,19 @@ struct ProjectivePipelineRoutines {
} // namespace
static void CompleteReconstructionLogProress(ProgressUpdateCallback *update_callback,
static void CompleteReconstructionLogProress(
ProgressUpdateCallback *update_callback,
double progress,
const char *step = NULL)
{
if(update_callback) {
const char *step = NULL) {
if (update_callback) {
char message[256];
if(step)
snprintf(message, sizeof(message), "Completing solution %d%% | %s", (int)(progress*100), step);
if (step)
snprintf(message, sizeof(message), "Completing solution %d%% | %s",
(int)(progress*100), step);
else
snprintf(message, sizeof(message), "Completing solution %d%%", (int)(progress*100));
snprintf(message, sizeof(message), "Completing solution %d%%",
(int)(progress*100));
update_callback->invoke(progress, message);
}
@@ -207,7 +213,8 @@ void InternalCompleteReconstruction(
if (reconstructed_markers.size() >= 5) {
CompleteReconstructionLogProress(update_callback,
(double)tot_resects/(max_image));
if (PipelineRoutines::Resect(options, reconstructed_markers, reconstruction, false)) {
if (PipelineRoutines::Resect(options, reconstructed_markers,
reconstruction, false)) {
num_resects++;
tot_resects++;
LG << "Ran Resect() for image " << image;
@@ -243,7 +250,8 @@ void InternalCompleteReconstruction(
if (reconstructed_markers.size() >= 5) {
CompleteReconstructionLogProress(update_callback,
(double)tot_resects/(max_image));
if (PipelineRoutines::Resect(options, reconstructed_markers, reconstruction, true)) {
if (PipelineRoutines::Resect(options, reconstructed_markers,
reconstruction, true)) {
num_resects++;
LG << "Ran final Resect() for image " << image;
} else {
@@ -260,9 +268,10 @@ void InternalCompleteReconstruction(
}
template<typename PipelineRoutines>
double InternalReprojectionError(const Tracks &image_tracks,
const typename PipelineRoutines::Reconstruction &reconstruction,
const CameraIntrinsics &intrinsics) {
double InternalReprojectionError(
const Tracks &image_tracks,
const typename PipelineRoutines::Reconstruction &reconstruction,
const CameraIntrinsics &intrinsics) {
int num_skipped = 0;
int num_reprojected = 0;
double total_error = 0.0;

9
extern/libmv/libmv/simple_pipeline/pipeline.h vendored
View File

@@ -49,10 +49,11 @@ namespace libmv {
\sa EuclideanResect, EuclideanIntersect, EuclideanBundle
*/
void EuclideanCompleteReconstruction(const ReconstructionOptions &options,
const Tracks &tracks,
EuclideanReconstruction *reconstruction,
ProgressUpdateCallback *update_callback = NULL);
void EuclideanCompleteReconstruction(
const ReconstructionOptions &options,
const Tracks &tracks,
EuclideanReconstruction *reconstruction,
ProgressUpdateCallback *update_callback = NULL);
/*!
Estimate camera matrices and homogeneous 3D coordinates for all frames and

14
extern/libmv/libmv/simple_pipeline/reconstruction.h vendored
View File

@@ -27,14 +27,14 @@
namespace libmv {
struct ReconstructionOptions {
// threshold value of reconstruction error which is still considered successful
// if reconstruction error bigger than this value, fallback reconstruction
// algorithm would be used (if enabled)
double success_threshold;
// threshold value of reconstruction error which is still considered successful
// if reconstruction error bigger than this value, fallback reconstruction
// algorithm would be used (if enabled)
double success_threshold;
// use fallback reconstruction algorithm in cases main reconstruction algorithm
// failed to reconstruct
bool use_fallback_reconstruction;
// use fallback reconstruction algorithm in cases main reconstruction algorithm
// failed to reconstruct
bool use_fallback_reconstruction;
};
/*!

23
extern/libmv/libmv/simple_pipeline/resect.cc vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2011 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -18,6 +18,8 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/simple_pipeline/resect.h"
#include <cstdio>
#include "libmv/base/vector.h"
@@ -27,7 +29,6 @@
#include "libmv/multiview/projection.h"
#include "libmv/numeric/numeric.h"
#include "libmv/numeric/levenberg_marquardt.h"
#include "libmv/simple_pipeline/resect.h"
#include "libmv/simple_pipeline/reconstruction.h"
#include "libmv/simple_pipeline/tracks.h"
@@ -108,13 +109,13 @@ bool EuclideanResect(const ReconstructionOptions &options,
double success_threshold = std::numeric_limits<double>::max();
if(options.use_fallback_reconstruction)
if (options.use_fallback_reconstruction)
success_threshold = options.success_threshold;
if (0 || !euclidean_resection::EuclideanResection(points_2d, points_3d, &R, &t,
euclidean_resection::RESECTION_EPNP,
success_threshold))
{
if (0 || !euclidean_resection::EuclideanResection(
points_2d, points_3d, &R, &t,
euclidean_resection::RESECTION_EPNP,
success_threshold)) {
// printf("Resection for image %d failed\n", markers[0].image);
LG << "Resection for image " << markers[0].image << " failed;"
<< " trying fallback projective resection.";
@@ -172,7 +173,7 @@ bool EuclideanResect(const ReconstructionOptions &options,
Solver solver(resect_cost);
Solver::SolverParameters params;
Solver::Results results = solver.minimize(params, &dRt);
/* Solver::Results results = */ solver.minimize(params, &dRt);
LG << "LM found incremental rotation: " << dRt.head<3>().transpose();
// TODO(keir): Check results to ensure clean termination.
@@ -264,7 +265,7 @@ bool ProjectiveResect(const vector<Marker> &markers,
Solver solver(resect_cost);
Solver::SolverParameters params;
Solver::Results results = solver.minimize(params, &vector_P);
/* Solver::Results results = */ solver.minimize(params, &vector_P);
// TODO(keir): Check results to ensure clean termination.
// Unpack the projection matrix.

6
extern/libmv/libmv/simple_pipeline/resect.h vendored
View File

@@ -1,15 +1,15 @@
// Copyright (c) 2011 libmv authors.
//
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to
// deal in the Software without restriction, including without limitation the
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
// sell copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

8
extern/libmv/libmv/simple_pipeline/tracks.cc vendored
View File

@@ -18,12 +18,13 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/simple_pipeline/tracks.h"
#include <algorithm>
#include <vector>
#include <iterator>
#include "libmv/numeric/numeric.h"
#include "libmv/simple_pipeline/tracks.h"
namespace libmv {
@@ -82,7 +83,8 @@ vector<Marker> Tracks::MarkersInBothImages(int image1, int image2) const {
return markers;
}
vector<Marker> Tracks::MarkersForTracksInBothImages(int image1, int image2) const {
vector<Marker> Tracks::MarkersForTracksInBothImages(int image1,
int image2) const {
std::vector<int> image1_tracks;
std::vector<int> image2_tracks;
@@ -106,7 +108,7 @@ vector<Marker> Tracks::MarkersForTracksInBothImages(int image1, int image2) cons
vector<Marker> markers;
for (int i = 0; i < markers_.size(); ++i) {
if ((markers_[i].image == image1 || markers_[i].image == image2) &&
std::binary_search(intersection.begin(),intersection.end(),
std::binary_search(intersection.begin(), intersection.end(),
markers_[i].track)) {
markers.push_back(markers_[i]);
}

19
extern/libmv/libmv/tracking/brute_region_tracker.cc vendored
View File

@@ -63,13 +63,13 @@ void *aligned_malloc(int size, int alignment) {
#elif __FreeBSD__
void *result;
if(posix_memalign(&result, alignment, size)) {
if (posix_memalign(&result, alignment, size)) {
// non-zero means allocation error
// either no allocation or bad alignment value
return NULL;
}
return result;
#else // This is for Linux.
#else // This is for Linux.
return memalign(alignment, size);
#endif
}
@@ -107,7 +107,7 @@ bool RegionIsInBounds(const FloatImage &image1,
#ifdef __SSE2__
// Compute the sub of absolute differences between the arrays "a" and "b".
// Compute the sub of absolute differences between the arrays "a" and "b".
// The array "a" is assumed to be 16-byte aligned, while "b" is not. The
// result is returned as the first and third elements of __m128i if
// interpreted as a 4-element 32-bit integer array. The SAD is the sum of the
@@ -122,7 +122,7 @@ inline static __m128i SumOfAbsoluteDifferencesContiguousSSE(
unsigned int size,
__m128i sad) {
// Do the bulk of the work as 16-way integer operations.
for(unsigned int j = 0; j < size / 16; j++) {
for (unsigned int j = 0; j < size / 16; j++) {
sad = _mm_add_epi32(sad, _mm_sad_epu8( _mm_load_si128 ((__m128i*)(a + 16 * j)),
_mm_loadu_si128((__m128i*)(b + 16 * j))));
}
@@ -307,7 +307,8 @@ bool BruteRegionTracker::Track(const FloatImage &image1,
// Convert the search area directly to bytes without sampling.
unsigned char *search_area;
int search_area_stride;
FloatArrayToByteArrayWithPadding(image_and_gradient2, &search_area, &search_area_stride);
FloatArrayToByteArrayWithPadding(image_and_gradient2, &search_area,
&search_area_stride);
// Try all possible locations inside the search area. Yes, everywhere.
int best_i = -1, best_j = -1, best_sad = INT_MAX;
@@ -361,9 +362,11 @@ bool BruteRegionTracker::Track(const FloatImage &image1,
// Compute the Pearson product-moment correlation coefficient to check
// for sanity.
double correlation = PearsonProductMomentCorrelation(image_and_gradient1_sampled,
image_and_gradient2_sampled,
pattern_width);
double correlation = PearsonProductMomentCorrelation(
image_and_gradient1_sampled,
image_and_gradient2_sampled,
pattern_width);
LG << "Final correlation: " << correlation;
if (correlation < minimum_correlation) {

2
extern/libmv/libmv/tracking/brute_region_tracker.h vendored
View File

@@ -30,7 +30,7 @@ struct BruteRegionTracker : public RegionTracker {
BruteRegionTracker()
: half_window_size(4),
minimum_correlation(0.78) {}
virtual ~BruteRegionTracker() {}
// Tracker interface.

18
extern/libmv/libmv/tracking/esm_region_tracker.cc vendored
View File

@@ -72,7 +72,7 @@ bool EsmRegionTracker::Track(const FloatImage &image1,
<< ", hw=" << half_window_size << ".";
return false;
}
// XXX
// TODO(keir): Delete the block between the XXX's once the planar tracker is
// integrated into blender.
@@ -101,7 +101,7 @@ bool EsmRegionTracker::Track(const FloatImage &image1,
options.max_iterations = 20;
options.sigma = sigma;
options.use_esm = true;
TrackRegionResult result;
TrackRegion(image1, image2, xx1, yy1, options, xx2, yy2, &result);
@@ -129,7 +129,7 @@ bool EsmRegionTracker::Track(const FloatImage &image1,
&image_and_gradient1_sampled);
// Step 0: Initialize delta = 0.01.
//
//
// Ignored for my "normal" LM loop.
// Step 1: Warp I with W(x, p) to compute I(W(x; p).
@@ -179,8 +179,7 @@ bool EsmRegionTracker::Track(const FloatImage &image1,
}
}
double tau = 1e-4, eps1, eps2, eps3;
eps1 = eps2 = eps3 = 1e-15;
double tau = 1e-4;
double mu = tau * std::max(H_image1(0, 0), H_image1(1, 1));
double nu = M_E;
@@ -201,7 +200,7 @@ bool EsmRegionTracker::Track(const FloatImage &image1,
image_and_gradient1_sampled(r, c, 2));
Vec2 g2(image_and_gradient2_sampled[current_image](r, c, 1),
image_and_gradient2_sampled[current_image](r, c, 2));
Vec2 g = g1 + g2; // Should be / 2.0, but do that outside the loop.
Vec2 g = g1 + g2; // Should be / 2.0, but do that outside the loop.
H += g * g.transpose();
}
}
@@ -282,9 +281,10 @@ bool EsmRegionTracker::Track(const FloatImage &image1,
if (d.squaredNorm() < min_update_squared_distance) {
// Compute the Pearson product-moment correlation coefficient to check
// for sanity.
double correlation = PearsonProductMomentCorrelation(image_and_gradient1_sampled,
image_and_gradient2_sampled[new_image],
width);
double correlation = PearsonProductMomentCorrelation(
image_and_gradient1_sampled,
image_and_gradient2_sampled[new_image],
width);
LG << "Final correlation: " << correlation;
// Note: Do the comparison here to handle nan's correctly (since all

2
extern/libmv/libmv/tracking/esm_region_tracker.h vendored
View File

@@ -40,7 +40,7 @@ struct EsmRegionTracker : public RegionTracker {
min_update_squared_distance(1e-4),
sigma(0.9),
minimum_correlation(0.78) {}
virtual ~EsmRegionTracker() {}
// Tracker interface.

2
extern/libmv/libmv/tracking/hybrid_region_tracker.h vendored
View File

@@ -34,7 +34,7 @@ class HybridRegionTracker : public RegionTracker {
RegionTracker *fine_tracker)
: coarse_tracker_(coarse_tracker),
fine_tracker_(fine_tracker) {}
virtual ~HybridRegionTracker() {}
// Tracker interface.

6
extern/libmv/libmv/tracking/klt_region_tracker.cc vendored
View File

@@ -18,8 +18,9 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/logging/logging.h"
#include "libmv/tracking/klt_region_tracker.h"
#include "libmv/logging/logging.h"
#include "libmv/image/image.h"
#include "libmv/image/convolve.h"
#include "libmv/image/sample.h"
@@ -144,7 +145,8 @@ bool KltRegionTracker::Track(const FloatImage &image1,
LG << "Determinant " << determinant << " is too small; failing tracking.";
return false;
}
LG << "x=" << *x2 << ", y=" << *y2 << ", dx=" << dx << ", dy=" << dy << ", det=" << determinant;
LG << "x=" << *x2 << ", y=" << *y2 << ", dx=" << dx << ", dy=" << dy
<< ", det=" << determinant;
// If the update is small, then we probably found the target.
if (dx * dx + dy * dy < min_update_squared_distance) {

2
extern/libmv/libmv/tracking/klt_region_tracker.h vendored
View File

@@ -33,7 +33,7 @@ struct KltRegionTracker : public RegionTracker {
min_determinant(1e-6),
min_update_squared_distance(1e-6),
sigma(0.9) {}
virtual ~KltRegionTracker() {}
// Tracker interface.

7
extern/libmv/libmv/tracking/lmicklt_region_tracker.cc vendored
View File

@@ -79,7 +79,7 @@ bool LmickltRegionTracker::Track(const FloatImage &image1,
<< ", hw=" << half_window_size << ".";
return false;
}
int width = 2 * half_window_size + 1;
// TODO(keir): Avoid recomputing gradients for e.g. the pyramid tracker.
@@ -100,7 +100,7 @@ bool LmickltRegionTracker::Track(const FloatImage &image1,
&image_and_gradient1_sampled);
// Step 0: Initialize delta = 0.01.
//
//
// Ignored for my "normal" LM loop.
double reasonable_error =
@@ -152,8 +152,7 @@ bool LmickltRegionTracker::Track(const FloatImage &image1,
}
}
double tau = 1e-3, eps1, eps2, eps3;
eps1 = eps2 = eps3 = 1e-15;
double tau = 1e-3;
double mu = tau * std::max(H(0, 0), H(1, 1));
double nu = 2.0;

2
extern/libmv/libmv/tracking/lmicklt_region_tracker.h vendored
View File

@@ -40,7 +40,7 @@ struct LmickltRegionTracker : public RegionTracker {
min_determinant(1e-6),
min_update_squared_distance(1e-6),
sigma(0.9) {}
virtual ~LmickltRegionTracker() {}
// Tracker interface.

3
extern/libmv/libmv/tracking/pyramid_region_tracker.cc vendored
View File

@@ -18,13 +18,14 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/tracking/pyramid_region_tracker.h"
#include <vector>
#include "libmv/image/convolve.h"
#include "libmv/image/image.h"
#include "libmv/image/sample.h"
#include "libmv/logging/logging.h"
#include "libmv/tracking/pyramid_region_tracker.h"
namespace libmv {

4
extern/libmv/libmv/tracking/retrack_region_tracker.cc vendored
View File

@@ -18,11 +18,11 @@
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
// IN THE SOFTWARE.
#include "libmv/tracking/retrack_region_tracker.h"
#include <cmath>
#include <vector>
#include "libmv/tracking/retrack_region_tracker.h"
namespace libmv {
bool RetrackRegionTracker::Track(const FloatImage &image1,

70
extern/libmv/libmv/tracking/track_region.cc vendored
View File

@@ -60,6 +60,8 @@ struct JetOps {
}
static void ScaleDerivative(double scale_by, T *value) {
// For double, there is no derivative to scale.
(void) scale_by; // Ignored.
(void) value; // Ignored.
}
};
@@ -85,7 +87,9 @@ struct Chain {
const FunctionType dfdx[kNumArgs],
const ArgumentType x[kNumArgs]) {
// In the default case of scalars, there's nothing to do since there are no
// derivatives to propagate.
// derivatives to propagate.
(void) dfdx; // Ignored.
(void) x; // Ignored.
return f;
}
};
@@ -276,8 +280,8 @@ class PixelDifferenceCostFunctor {
int num_samples_y,
const Warp &warp)
: options_(options),
image_and_gradient1_(image_and_gradient1),
image_and_gradient2_(image_and_gradient2),
image_and_gradient1_(image_and_gradient1),
image_and_gradient2_(image_and_gradient2),
canonical_to_image1_(canonical_to_image1),
num_samples_x_(num_samples_x),
num_samples_y_(num_samples_y),
@@ -348,7 +352,7 @@ class PixelDifferenceCostFunctor {
//
// Note that partial masks are not short circuited. To see why short
// circuiting produces bitwise-exact same results, consider that the
// residual for each pixel is
// residual for each pixel is
//
// residual = mask * (src - dst) ,
//
@@ -431,11 +435,10 @@ class PixelDifferenceCostFunctor {
return true;
}
// For normalized matching, the average and
// For normalized matching, the average and
template<typename T>
void ComputeNormalizingCoefficient(const T *warp_parameters,
T *dst_mean) const {
*dst_mean = T(0.0);
double num_samples = 0.0;
for (int r = 0; r < num_samples_y_; ++r) {
@@ -443,7 +446,7 @@ class PixelDifferenceCostFunctor {
// Use the pre-computed image1 position.
Vec2 image1_position(pattern_positions_(r, c, 0),
pattern_positions_(r, c, 1));
// Sample the mask early; if it's zero, this pixel has no effect. This
// allows early bailout from the expensive sampling that happens below.
double mask_value = 1.0;
@@ -484,28 +487,28 @@ class PixelDifferenceCostFunctor {
LG << "Normalization for dst:" << *dst_mean;
}
// TODO(keir): Consider also computing the cost here.
double PearsonProductMomentCorrelationCoefficient(
const double *warp_parameters) const {
for (int i = 0; i < Warp::NUM_PARAMETERS; ++i) {
VLOG(2) << "Correlation warp_parameters[" << i << "]: "
<< warp_parameters[i];
}
// TODO(keir): Consider also computing the cost here.
double PearsonProductMomentCorrelationCoefficient(
const double *warp_parameters) const {
for (int i = 0; i < Warp::NUM_PARAMETERS; ++i) {
VLOG(2) << "Correlation warp_parameters[" << i << "]: "
<< warp_parameters[i];
}
// The single-pass PMCC computation is somewhat numerically unstable, but
// it's sufficient for the tracker.
double sX = 0, sY = 0, sXX = 0, sYY = 0, sXY = 0;
// The single-pass PMCC computation is somewhat numerically unstable, but
// it's sufficient for the tracker.
double sX = 0, sY = 0, sXX = 0, sYY = 0, sXY = 0;
// Due to masking, it's important to account for fractional samples.
// For example, samples with a 50% mask are counted as a half sample.
double num_samples = 0;
// Due to masking, it's important to account for fractional samples.
// For example, samples with a 50% mask are counted as a half sample.
double num_samples = 0;
for (int r = 0; r < num_samples_y_; ++r) {
for (int c = 0; c < num_samples_x_; ++c) {
for (int r = 0; r < num_samples_y_; ++r) {
for (int c = 0; c < num_samples_x_; ++c) {
// Use the pre-computed image1 position.
Vec2 image1_position(pattern_positions_(r, c, 0),
pattern_positions_(r, c, 1));
double mask_value = 1.0;
if (options_.image1_mask != NULL) {
mask_value = pattern_mask_(r, c);
@@ -795,7 +798,7 @@ struct TranslationRotationWarp {
// Obtain the rotation via orthorgonal procrustes.
Mat2 correlation_matrix;
for (int i = 0; i < 4; ++i) {
correlation_matrix += q1.CornerRelativeToCentroid(i) *
correlation_matrix += q1.CornerRelativeToCentroid(i) *
q2.CornerRelativeToCentroid(i).transpose();
}
Mat2 R = OrthogonalProcrustes(correlation_matrix);
@@ -863,7 +866,7 @@ struct TranslationRotationScaleWarp {
// Obtain the rotation via orthorgonal procrustes.
Mat2 correlation_matrix;
for (int i = 0; i < 4; ++i) {
correlation_matrix += q1.CornerRelativeToCentroid(i) *
correlation_matrix += q1.CornerRelativeToCentroid(i) *
q2.CornerRelativeToCentroid(i).transpose();
}
Mat2 R = OrthogonalProcrustes(correlation_matrix);
@@ -938,7 +941,7 @@ struct AffineWarp {
Vec2 v1 = q1.CornerRelativeToCentroid(i);
Vec2 v2 = q2.CornerRelativeToCentroid(i);
Q1.row(2 * i + 0) << v1[0], v1[1], 0, 0 ;
Q1.row(2 * i + 0) << v1[0], v1[1], 0, 0;
Q1.row(2 * i + 1) << 0, 0, v1[0], v1[1];
Q2(2 * i + 0) = v2[0];
@@ -1039,6 +1042,9 @@ struct HomographyWarp {
void PickSampling(const double *x1, const double *y1,
const double *x2, const double *y2,
int *num_samples_x, int *num_samples_y) {
(void) x2; // Ignored.
(void) y2; // Ignored.
Vec2 a0(x1[0], y1[0]);
Vec2 a1(x1[1], y1[1]);
Vec2 a2(x1[2], y1[2]);
@@ -1074,6 +1080,10 @@ void PickSampling(const double *x1, const double *y1,
bool SearchAreaTooBigForDescent(const FloatImage &image2,
const double *x2, const double *y2) {
// TODO(keir): Check the bounds and enable only when it makes sense.
(void) image2; // Ignored.
(void) x2; // Ignored.
(void) y2; // Ignored.
return true;
}
@@ -1150,12 +1160,12 @@ void CreateBrutePattern(const double *x1, const double *y1,
inverse_warp.Forward(inverse_warp.parameters,
dst_x, dst_y,
&src_x, &src_y);
if (PointInQuad(x1, y1, src_x, src_y)) {
(*pattern)(i, j) = SampleLinear(image1, src_y, src_x);
(*mask)(i, j) = 1.0;
if (image1_mask) {
(*mask)(i, j) = SampleLinear(*image1_mask, src_y, src_x);;
(*mask)(i, j) = SampleLinear(*image1_mask, src_y, src_x);
}
} else {
(*pattern)(i, j) = 0.0;
@@ -1218,7 +1228,7 @@ bool BruteTranslationOnlyInitialize(const FloatImage &image1,
//
// TODO(keir): There are a number of possible optimizations here. One choice
// is to make a grid and only try one out of every N possible samples.
//
//
// Another, slightly more clever idea, is to compute some sort of spatial
// frequency distribution of the pattern patch. If the spatial resolution is
// high (e.g. a grating pattern or fine lines) then checking every possible
@@ -1375,7 +1385,7 @@ void TemplatedTrackRegion(const FloatImage &image1,
num_samples_x,
num_samples_y,
warp);
problem.AddResidualBlock(
problem.AddResidualBlock(
new ceres::AutoDiffCostFunction<
PixelDifferenceCostFunctor<Warp>,
ceres::DYNAMIC,

3
extern/libmv/libmv/tracking/trklt_region_tracker.cc vendored
View File

@@ -158,7 +158,8 @@ bool TrkltRegionTracker::Track(const FloatImage &image1,
LG << "Determinant " << determinant << " is too small; failing tracking.";
return false;
}
LG << "x=" << *x2 << ", y=" << *y2 << ", dx=" << d[0] << ", dy=" << d[1] << ", det=" << determinant;
LG << "x=" << *x2 << ", y=" << *y2 << ", dx=" << d[0] << ", dy=" << d[1]
<< ", det=" << determinant;
// If the update is small, then we probably found the target.

3
extern/libmv/third_party/fast/fast_10.c vendored
View File

@@ -1,8 +1,7 @@
/*This is mechanically generated code*/
#include <stdlib.h>
typedef struct { int x, y; } xy;
typedef unsigned char byte;
#include "fast.h"
int fast10_corner_score(const byte* p, const int pixel[], int bstart)
{

3
extern/libmv/third_party/fast/fast_11.c vendored
View File

@@ -1,8 +1,7 @@
/*This is mechanically generated code*/
#include <stdlib.h>
typedef struct { int x, y; } xy;
typedef unsigned char byte;
#include "fast.h"
int fast11_corner_score(const byte* p, const int pixel[], int bstart)
{

3
extern/libmv/third_party/fast/fast_12.c vendored
View File

@@ -1,8 +1,7 @@
/*This is mechanically generated code*/
#include <stdlib.h>
typedef struct { int x, y; } xy;
typedef unsigned char byte;
#include "fast.h"
int fast12_corner_score(const byte* p, const int pixel[], int bstart)
{

3
extern/libmv/third_party/fast/fast_9.c vendored
View File

@@ -1,8 +1,7 @@
/*This is mechanically generated code*/
#include <stdlib.h>
typedef struct { int x, y; } xy;
typedef unsigned char byte;
#include "fast.h"
int fast9_corner_score(const byte* p, const int pixel[], int bstart)
{

1
extern/libmv/third_party/gflags/gflags_completions.cc vendored
View File

@@ -57,6 +57,7 @@
#include <utility>
#include <vector>
#include "gflags/gflags_completions.h"
#include "gflags/gflags.h"
#include "util.h"

4
extern/libmv/third_party/glog/src/logging.cc vendored
View File

@@ -1464,6 +1464,10 @@ void LogToStderr() {
namespace base {
namespace internal {
/* Put prototypes here to suppress strict compiler warnings */
bool GetExitOnDFatal();
void SetExitOnDFatal(bool value);
bool GetExitOnDFatal() {
MutexLock l(&log_mutex);
return exit_on_dfatal;

2
extern/libmv/third_party/glog/src/utilities.cc vendored
View File

@@ -84,7 +84,7 @@ static void DebugWriteToStderr(const char* data, void *) {
}
}
void DebugWriteToString(const char* data, void *arg) {
static void DebugWriteToString(const char* data, void *arg) {
reinterpret_cast<string*>(arg)->append(data);
}

6
extern/libmv/third_party/glog/src/vlog_is_on.cc vendored
View File

@@ -62,6 +62,12 @@ _START_GOOGLE_NAMESPACE_
namespace glog_internal_namespace_ {
// Put protytype here to suppress strict compiler flags
GOOGLE_GLOG_DLL_DECL bool SafeFNMatch_(const char* pattern,
size_t patt_len,
const char* str,
size_t str_len);
// Implementation of fnmatch that does not need 0-termination
// of arguments and does not allocate any memory,
// but we only support "*" and "?" wildcards, not the "[...]" patterns.