diff --git a/.cproject b/.cproject
index d5a6ca4d4..62330cdbc 100644
--- a/.cproject
+++ b/.cproject
@@ -518,6 +518,22 @@
true
true
+
+ make
+ -j5
+ testInvDepthCamera3.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testTriangulation.run
+ true
+ true
+ true
+
make
-j2
@@ -662,62 +678,6 @@
false
true
-
- make
- -j5
- testPoseRTV.run
- true
- true
- true
-
-
- make
- -j5
- testIMUSystem.run
- true
- true
- true
-
-
- make
- -j5
- testBTree.run
- true
- true
- true
-
-
- make
- -j5
- testDSF.run
- true
- true
- true
-
-
- make
- -j5
- testDSFVector.run
- true
- true
- true
-
-
- make
- -j5
- testFixedVector.run
- true
- true
- true
-
-
- make
- -j5
- testDSFMap.run
- true
- true
- true
-
make
-j2
@@ -758,50 +718,74 @@
true
true
-
+
make
-j5
- schedulingExample.run
+ testAHRS.run
true
true
true
-
+
make
-j5
- testCSP.run
+ testImuFactor.run
true
true
true
-
+
make
-j5
- testScheduler.run
+ testInvDepthFactor3.run
true
true
true
-
+
make
-j5
- schedulingQuals12.run
+ testMultiProjectionFactor.run
true
true
true
-
+
make
-j5
- testSudoku.run
+ testPoseRotationPrior.run
true
true
true
-
+
make
-j5
- schedulingQuals13.run
+ testPoseTranslationPrior.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testReferenceFrameFactor.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testSmartProjectionFactor.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testTSAMFactors.run
true
true
true
@@ -902,78 +886,22 @@
true
true
-
+
make
-j5
- testInvDepthFactor3.run
+ testCombinedImuFactor.run
true
true
true
-
+
make
-j5
- testPoseTranslationPrior.run
- true
- true
- true
-
-
- make
- -j5
- testPoseRotationPrior.run
- true
- true
- true
-
-
- make
- -j5
- testReferenceFrameFactor.run
- true
- true
- true
-
-
- make
- -j5
- testAHRS.run
- true
- true
- true
-
-
- make
- -j8
testImuFactor.run
true
true
true
-
- make
- -j5
- testMultiProjectionFactor.run
- true
- true
- true
-
-
- make
- -j5
- testSmartProjectionFactor.run
- true
- true
- true
-
-
- make
- -j5
- testTSAMFactors.run
- true
- true
- true
-
make
-j5
@@ -1126,14 +1054,6 @@
true
true
-
- make
- -j5
- testParticleFactor.run
- true
- true
- true
-
make
-j2
@@ -1318,18 +1238,42 @@
true
true
-
+
make
-j5
- testImuFactor.run
+ testBTree.run
true
true
true
-
+
make
-j5
- testCombinedImuFactor.run
+ testDSF.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testDSFMap.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testDSFVector.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testFixedVector.run
true
true
true
@@ -2081,26 +2025,10 @@
true
true
-
+
make
-j5
- timeCalibratedCamera.run
- true
- true
- true
-
-
- make
- -j5
- timePinholeCamera.run
- true
- true
- true
-
-
- make
- -j5
- timeStereoCamera.run
+ testCal3Unified.run
true
true
true
@@ -2185,6 +2113,62 @@
true
true
+
+ make
+ -j5
+ testIMUSystem.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testPoseRTV.run
+ true
+ true
+ true
+
+
+ make
+ -j1
+ testDiscreteBayesTree.run
+ true
+ false
+ true
+
+
+ make
+ -j5
+ testDiscreteConditional.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testDiscreteFactor.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testDiscreteFactorGraph.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testDiscreteMarginals.run
+ true
+ true
+ true
+
make
-j2
@@ -2225,42 +2209,50 @@
true
true
-
+
make
-j5
- testDiscreteFactor.run
+ schedulingExample.run
true
true
true
-
- make
- -j1
- testDiscreteBayesTree.run
- true
- false
- true
-
-
+
make
-j5
- testDiscreteFactorGraph.run
+ schedulingQuals12.run
true
true
true
-
+
make
-j5
- testDiscreteConditional.run
+ schedulingQuals13.run
true
true
true
-
+
make
-j5
- testDiscreteMarginals.run
+ testCSP.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testScheduler.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ testSudoku.run
true
true
true
@@ -2281,22 +2273,6 @@
true
true
-
- make
- -j5
- testInvDepthCamera3.run
- true
- true
- true
-
-
- make
- -j5
- testTriangulation.run
- true
- true
- true
-
make
-j5
@@ -2529,6 +2505,14 @@
true
true
+
+ make
+ -j5
+ testParticleFactor.run
+ true
+ true
+ true
+
make
-j2
@@ -2817,6 +2801,14 @@
true
true
+
+ make
+ -j5
+ SFMExample_SmartFactor.run
+ true
+ true
+ true
+
make
-j5
@@ -2857,7 +2849,39 @@
true
true
-
+
+ make
+ -j4
+ testImuFactor.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ timeCalibratedCamera.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ timePinholeCamera.run
+ true
+ true
+ true
+
+
+ make
+ -j5
+ timeStereoCamera.run
+ true
+ true
+ true
+
+
make
-j5
timeLago.run
@@ -2865,10 +2889,10 @@
true
true
-
+
make
- -j4
- testImuFactor.run
+ -j5
+ timePose3.run
true
true
true
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 933f2083e..004ded5e4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -34,7 +34,7 @@ if(${CMAKE_SOURCE_DIR} STREQUAL ${CMAKE_BINARY_DIR})
message(FATAL_ERROR "In-source builds not allowed. Please make a new directory (called a build directory) and run CMake from there. You may need to remove CMakeCache.txt. ")
endif()
-# See whether gtsam_unstable is available (it will be present only if we're using an SVN checkout)
+# See whether gtsam_unstable is available (it will be present only if we're using a git checkout)
if(EXISTS "${PROJECT_SOURCE_DIR}/gtsam_unstable" AND IS_DIRECTORY "${PROJECT_SOURCE_DIR}/gtsam_unstable")
set(GTSAM_UNSTABLE_AVAILABLE 1)
else()
@@ -46,7 +46,6 @@ endif()
# Set up options
# Configurable Options
-option(GTSAM_BUILD_TIMING "Enable/Disable building of timing scripts" OFF) # These do not currently work
if(GTSAM_UNSTABLE_AVAILABLE)
option(GTSAM_BUILD_UNSTABLE "Enable/Disable libgtsam_unstable" ON)
endif()
@@ -94,7 +93,7 @@ if(MSVC)
# Disable autolinking
if(NOT Boost_USE_STATIC_LIBS)
add_definitions(-DBOOST_ALL_NO_LIB)
- add_definitions(-DBOOST_ALL_DYN_LINK)
+ add_definitions(-DBOOST_ALL_DYN_LINK)
endif()
endif()
@@ -156,21 +155,21 @@ find_package(GooglePerfTools)
###############################################################################
# Find MKL
-if(GTSAM_USE_EIGEN_MKL)
- find_package(MKL)
+find_package(MKL)
- if(MKL_FOUND AND GTSAM_WITH_EIGEN_MKL)
- set(GTSAM_USE_EIGEN_MKL 1) # This will go into config.h
- set(EIGEN_USE_MKL_ALL 1) # This will go into config.h - it makes Eigen use MKL
- include_directories(${MKL_INCLUDE_DIR})
- list(APPEND GTSAM_ADDITIONAL_LIBRARIES ${MKL_LIBRARIES})
- endif()
+if(MKL_FOUND AND GTSAM_WITH_EIGEN_MKL)
+ set(GTSAM_USE_EIGEN_MKL 1) # This will go into config.h
+ set(EIGEN_USE_MKL_ALL 1) # This will go into config.h - it makes Eigen use MKL
+ include_directories(${MKL_INCLUDE_DIR})
+ list(APPEND GTSAM_ADDITIONAL_LIBRARIES ${MKL_LIBRARIES})
+else()
+ set(GTSAM_USE_EIGEN_MKL 0)
+ set(EIGEN_USE_MKL_ALL 0)
endif()
-
###############################################################################
# Find OpenMP (if we're also using MKL)
-if(GTSAM_USE_EIGEN_MKL AND GTSAM_USE_EIGEN_MKL_OPENMP AND GTSAM_USE_EIGEN_MKL)
+if(GTSAM_WITH_EIGEN_MKL AND GTSAM_USE_EIGEN_MKL_OPENMP AND GTSAM_USE_EIGEN_MKL)
find_package(OpenMP)
if(OPENMP_FOUND AND GTSAM_USE_EIGEN_MKL AND GTSAM_WITH_EIGEN_MKL_OPENMP)
@@ -304,6 +303,9 @@ add_subdirectory(tests)
# Build examples
add_subdirectory(examples)
+# Build timing
+add_subdirectory(timing)
+
# Matlab toolbox
if (GTSAM_INSTALL_MATLAB_TOOLBOX)
add_subdirectory(matlab)
@@ -361,7 +363,7 @@ message(STATUS "================ Configuration Options ======================"
message(STATUS "Build flags ")
print_config_flag(${GTSAM_BUILD_TESTS} "Build Tests ")
print_config_flag(${GTSAM_BUILD_EXAMPLES_ALWAYS} "Build examples with 'make all' ")
-print_config_flag(${GTSAM_BUILD_TIMING} "Build Timing scripts ")
+print_config_flag(${GTSAM_BUILD_TIMING_ALWAYS} "Build timing scripts with 'make all'")
if (DOXYGEN_FOUND)
print_config_flag(${GTSAM_BUILD_DOCS} "Build Docs ")
endif()
diff --git a/cmake/GTSAMCMakeToolsConfig.cmake b/cmake/GTSAMCMakeToolsConfig.cmake
index c79a2f5c2..126bcd2d9 100644
--- a/cmake/GTSAMCMakeToolsConfig.cmake
+++ b/cmake/GTSAMCMakeToolsConfig.cmake
@@ -1,3 +1,4 @@
# This config file modifies CMAKE_MODULE_PATH so that the GTSAM-CMakeTools files may be included
-
+
+set(GTSAM_CMAKE_TOOLS_DIR "${CMAKE_CURRENT_LIST_DIR}")
list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}")
diff --git a/cmake/GtsamMatlabWrap.cmake b/cmake/GtsamMatlabWrap.cmake
index 5d018376b..ba616b025 100644
--- a/cmake/GtsamMatlabWrap.cmake
+++ b/cmake/GtsamMatlabWrap.cmake
@@ -28,13 +28,13 @@ endif()
# finding the LaTeX mex program (totally unrelated to MATLAB Mex) when LaTeX is
# on the system path.
list(REVERSE matlab_bin_directories) # Reverse list so the highest version (sorted alphabetically) is preferred
-find_program(mex_command ${mex_program_name}
+find_program(MEX_COMMAND ${mex_program_name}
PATHS ${matlab_bin_directories} ENV PATH
NO_DEFAULT_PATH)
-mark_as_advanced(FORCE mex_command)
+mark_as_advanced(FORCE MEX_COMMAND)
# Now that we have mex, trace back to find the Matlab installation root
-get_filename_component(mex_command "${mex_command}" REALPATH)
-get_filename_component(mex_path "${mex_command}" PATH)
+get_filename_component(MEX_COMMAND "${MEX_COMMAND}" REALPATH)
+get_filename_component(mex_path "${MEX_COMMAND}" PATH)
get_filename_component(MATLAB_ROOT "${mex_path}/.." ABSOLUTE)
set(MATLAB_ROOT "${MATLAB_ROOT}" CACHE PATH "Path to MATLAB installation root (e.g. /usr/local/MATLAB/R2012a)")
diff --git a/cmake/GtsamTesting.cmake b/cmake/GtsamTesting.cmake
index 19f487256..2e505c11e 100644
--- a/cmake/GtsamTesting.cmake
+++ b/cmake/GtsamTesting.cmake
@@ -38,21 +38,46 @@ endmacro()
#
# Add scripts that will serve as examples of how to use the library. A list of files or
# glob patterns is specified, and one executable will be created for each matching .cpp
-# file. These executables will not be installed. They are build with 'make all' if
+# file. These executables will not be installed. They are built with 'make all' if
# GTSAM_BUILD_EXAMPLES_ALWAYS is enabled. They may also be built with 'make examples'.
#
# Usage example:
# gtsamAddExamplesGlob("*.cpp" "BrokenExample.cpp" "gtsam;GeographicLib")
#
# Arguments:
-# globPatterns: The list of files or glob patterns from which to create unit tests, with
-# one test created for each cpp file. e.g. "*.cpp", or
+# globPatterns: The list of files or glob patterns from which to create examples, with
+# one program created for each cpp file. e.g. "*.cpp", or
# "A*.cpp;B*.cpp;MyExample.cpp".
# excludedFiles: A list of files or globs to exclude, e.g. "C*.cpp;BrokenExample.cpp". Pass
# an empty string "" if nothing needs to be excluded.
# linkLibraries: The list of libraries to link to.
macro(gtsamAddExamplesGlob globPatterns excludedFiles linkLibraries)
- gtsamAddExamplesGlob_impl("${globPatterns}" "${excludedFiles}" "${linkLibraries}")
+ gtsamAddExesGlob_impl("${globPatterns}" "${excludedFiles}" "${linkLibraries}" "examples" ${GTSAM_BUILD_EXAMPLES_ALWAYS})
+endmacro()
+
+
+###############################################################################
+# Macro:
+#
+# gtsamAddTimingGlob(globPatterns excludedFiles linkLibraries)
+#
+# Add scripts that time aspects of the library. A list of files or
+# glob patterns is specified, and one executable will be created for each matching .cpp
+# file. These executables will not be installed. They are not built with 'make all',
+# but they may be built with 'make timing'.
+#
+# Usage example:
+# gtsamAddTimingGlob("*.cpp" "DisabledTimingScript.cpp" "gtsam;GeographicLib")
+#
+# Arguments:
+# globPatterns: The list of files or glob patterns from which to create programs, with
+# one program created for each cpp file. e.g. "*.cpp", or
+# "A*.cpp;B*.cpp;MyExample.cpp".
+# excludedFiles: A list of files or globs to exclude, e.g. "C*.cpp;BrokenExample.cpp". Pass
+# an empty string "" if nothing needs to be excluded.
+# linkLibraries: The list of libraries to link to.
+macro(gtsamAddTimingGlob globPatterns excludedFiles linkLibraries)
+ gtsamAddExesGlob_impl("${globPatterns}" "${excludedFiles}" "${linkLibraries}" "timing" ${GTSAM_BUILD_TIMING_ALWAYS})
endmacro()
@@ -63,6 +88,7 @@ enable_testing()
option(GTSAM_BUILD_TESTS "Enable/Disable building of tests" ON)
option(GTSAM_BUILD_EXAMPLES_ALWAYS "Build examples with 'make all' (build with 'make examples' if not)" ON)
+option(GTSAM_BUILD_TIMING_ALWAYS "Build timing scripts with 'make all' (build with 'make timing' if not" OFF)
# Add option for combining unit tests
if(MSVC OR XCODE_VERSION)
@@ -80,6 +106,9 @@ endif()
# Add examples target
add_custom_target(examples)
+# Add timing target
+add_custom_target(timing)
+
# Include obsolete macros - will be removed in the near future
include(GtsamTestingObsolete)
@@ -180,7 +209,7 @@ macro(gtsamAddTestsGlob_impl groupName globPatterns excludedFiles linkLibraries)
endmacro()
-macro(gtsamAddExamplesGlob_impl globPatterns excludedFiles linkLibraries)
+macro(gtsamAddExesGlob_impl globPatterns excludedFiles linkLibraries groupName buildWithAll)
# Get all script files
file(GLOB script_files ${globPatterns})
@@ -220,20 +249,22 @@ macro(gtsamAddExamplesGlob_impl globPatterns excludedFiles linkLibraries)
target_link_libraries(${script_name} ${linkLibraries})
# Add target dependencies
- add_dependencies(examples ${script_name})
+ add_dependencies(${groupName} ${script_name})
if(NOT MSVC AND NOT XCODE_VERSION)
add_custom_target(${script_name}.run ${EXECUTABLE_OUTPUT_PATH}${script_name})
endif()
# Add TOPSRCDIR
set_property(SOURCE ${script_src} APPEND PROPERTY COMPILE_DEFINITIONS "TOPSRCDIR=\"${PROJECT_SOURCE_DIR}\"")
-
- if(NOT GTSAM_BUILD_EXAMPLES_ALWAYS)
+
+ # Exclude from all or not - note weird variable assignment because we're in a macro
+ set(buildWithAll_on ${buildWithAll})
+ if(NOT buildWithAll_on)
# Exclude from 'make all' and 'make install'
- set_target_properties(${target_name} PROPERTIES EXCLUDE_FROM_ALL ON)
+ set_target_properties("${script_name}" PROPERTIES EXCLUDE_FROM_ALL ON)
endif()
# Configure target folder (for MSVC and Xcode)
- set_property(TARGET ${script_name} PROPERTY FOLDER "Examples")
+ set_property(TARGET ${script_name} PROPERTY FOLDER "${groupName}")
endforeach()
endmacro()
diff --git a/examples/Pose2SLAMExample_g2o.cpp b/examples/Pose2SLAMExample_g2o.cpp
index 3a4ee9cff..b46a53198 100644
--- a/examples/Pose2SLAMExample_g2o.cpp
+++ b/examples/Pose2SLAMExample_g2o.cpp
@@ -31,7 +31,7 @@ int main(const int argc, const char *argv[]) {
// Read graph from file
string g2oFile;
if (argc < 2)
- g2oFile = "../../examples/Data/noisyToyGraph.txt";
+ g2oFile = findExampleDataFile("noisyToyGraph.txt");
else
g2oFile = argv[1];
diff --git a/examples/Pose2SLAMExample_graph.cpp b/examples/Pose2SLAMExample_graph.cpp
index 449144fef..86be74c75 100644
--- a/examples/Pose2SLAMExample_graph.cpp
+++ b/examples/Pose2SLAMExample_graph.cpp
@@ -32,7 +32,8 @@ int main (int argc, char** argv) {
NonlinearFactorGraph::shared_ptr graph;
Values::shared_ptr initial;
SharedDiagonal model = noiseModel::Diagonal::Sigmas((Vector(3) << 0.05, 0.05, 5.0 * M_PI / 180.0));
- boost::tie(graph, initial) = load2D("../../examples/Data/w100.graph", model);
+ string graph_file = findExampleDataFile("w100.graph");
+ boost::tie(graph, initial) = load2D(graph_file, model);
initial->print("Initial estimate:\n");
// Add a Gaussian prior on first poses
diff --git a/examples/Pose2SLAMExample_lago.cpp b/examples/Pose2SLAMExample_lago.cpp
index 1f5d80d7b..bb434f3ce 100644
--- a/examples/Pose2SLAMExample_lago.cpp
+++ b/examples/Pose2SLAMExample_lago.cpp
@@ -32,7 +32,7 @@ int main(const int argc, const char *argv[]) {
// Read graph from file
string g2oFile;
if (argc < 2)
- g2oFile = "../../examples/Data/noisyToyGraph.txt";
+ g2oFile = findExampleDataFile("noisyToyGraph.txt");
else
g2oFile = argv[1];
diff --git a/examples/Pose2SLAMwSPCG.cpp b/examples/Pose2SLAMwSPCG.cpp
index 2b2f63064..8e8e28570 100644
--- a/examples/Pose2SLAMwSPCG.cpp
+++ b/examples/Pose2SLAMwSPCG.cpp
@@ -104,7 +104,7 @@ int main(int argc, char** argv) {
LevenbergMarquardtParams parameters;
parameters.verbosity = NonlinearOptimizerParams::ERROR;
parameters.verbosityLM = LevenbergMarquardtParams::LAMBDA;
- parameters.linearSolverType = NonlinearOptimizerParams::CONJUGATE_GRADIENT;
+ parameters.linearSolverType = NonlinearOptimizerParams::Iterative;
{
parameters.iterativeParams = boost::make_shared();
diff --git a/examples/RangeISAMExample_plaza2.cpp b/examples/RangeISAMExample_plaza2.cpp
index 29b02a271..252372f39 100644
--- a/examples/RangeISAMExample_plaza2.cpp
+++ b/examples/RangeISAMExample_plaza2.cpp
@@ -40,6 +40,7 @@
#include
#include
#include
+#include
// Standard headers, added last, so we know headers above work on their own
#include
@@ -59,9 +60,8 @@ namespace NM = gtsam::noiseModel;
typedef pair TimedOdometry;
list readOdometry() {
list odometryList;
- ifstream is("../../examples/Data/Plaza2_DR.txt");
- if (!is)
- throw runtime_error("../../examples/Data/Plaza2_DR.txt file not found");
+ string data_file = findExampleDataFile("Plaza2_DR.txt");
+ ifstream is(data_file.c_str());
while (is) {
double t, distance_traveled, delta_heading;
@@ -78,9 +78,8 @@ list readOdometry() {
typedef boost::tuple RangeTriple;
vector readTriples() {
vector triples;
- ifstream is("../../examples/Data/Plaza2_TD.txt");
- if (!is)
- throw runtime_error("../../examples/Data/Plaza2_TD.txt file not found");
+ string data_file = findExampleDataFile("Plaza2_TD.txt");
+ ifstream is(data_file.c_str());
while (is) {
double t, sender, receiver, range;
diff --git a/examples/SFMExample_SmartFactor.cpp b/examples/SFMExample_SmartFactor.cpp
index 904919d42..b3c5ee5fe 100644
--- a/examples/SFMExample_SmartFactor.cpp
+++ b/examples/SFMExample_SmartFactor.cpp
@@ -14,6 +14,7 @@
* @brief A structure-from-motion problem on a simulated dataset, using smart projection factor
* @author Duy-Nguyen Ta
* @author Jing Dong
+ * @author Frank Dellaert
*/
/**
@@ -28,11 +29,6 @@
// Camera observations of landmarks (i.e. pixel coordinates) will be stored as Point2 (x, y).
#include
-// Each variable in the system (poses and landmarks) must be identified with a unique key.
-// We can either use simple integer keys (1, 2, 3, ...) or symbols (X1, X2, L1).
-// Here we will use Symbols
-#include
-
// In GTSAM, measurement functions are represented as 'factors'.
// The factor we used here is SmartProjectionPoseFactor. Every smart factor represent a single landmark,
// The SmartProjectionPoseFactor only optimize the pose of camera, not the calibration,
@@ -65,8 +61,8 @@ using namespace std;
using namespace gtsam;
// Make the typename short so it looks much cleaner
-typedef gtsam::SmartProjectionPoseFactor
- SmartFactor;
+typedef gtsam::SmartProjectionPoseFactor SmartFactor;
/* ************************************************************************* */
int main(int argc, char* argv[]) {
@@ -75,93 +71,87 @@ int main(int argc, char* argv[]) {
Cal3_S2::shared_ptr K(new Cal3_S2(50.0, 50.0, 0.0, 50.0, 50.0));
// Define the camera observation noise model
- noiseModel::Isotropic::shared_ptr measurementNoise = noiseModel::Isotropic::Sigma(2, 1.0); // one pixel in u and v
+ noiseModel::Isotropic::shared_ptr measurementNoise =
+ noiseModel::Isotropic::Sigma(2, 1.0); // one pixel in u and v
- // Create the set of ground-truth landmarks
+ // Create the set of ground-truth landmarks and poses
vector points = createPoints();
-
- // Create the set of ground-truth poses
vector poses = createPoses();
// Create a factor graph
NonlinearFactorGraph graph;
- // A vector saved all Smart factors (for get landmark position after optimization)
- vector smartfactors_ptr;
-
// Simulated measurements from each camera pose, adding them to the factor graph
- for (size_t i = 0; i < points.size(); ++i) {
+ for (size_t j = 0; j < points.size(); ++j) {
// every landmark represent a single landmark, we use shared pointer to init the factor, and then insert measurements.
SmartFactor::shared_ptr smartfactor(new SmartFactor());
- for (size_t j = 0; j < poses.size(); ++j) {
+ for (size_t i = 0; i < poses.size(); ++i) {
// generate the 2D measurement
- SimpleCamera camera(poses[j], *K);
- Point2 measurement = camera.project(points[i]);
+ SimpleCamera camera(poses[i], *K);
+ Point2 measurement = camera.project(points[j]);
- // call add() function to add measurment into a single factor, here we need to add:
+ // call add() function to add measurement into a single factor, here we need to add:
// 1. the 2D measurement
// 2. the corresponding camera's key
// 3. camera noise model
// 4. camera calibration
- smartfactor->add(measurement, Symbol('x', j), measurementNoise, K);
+ smartfactor->add(measurement, i, measurementNoise, K);
}
- // save smartfactors to get landmark position
- smartfactors_ptr.push_back(smartfactor);
-
// insert the smart factor in the graph
graph.push_back(smartfactor);
}
// Add a prior on pose x0. This indirectly specifies where the origin is.
- noiseModel::Diagonal::shared_ptr poseNoise = noiseModel::Diagonal::Sigmas((Vector(6) << Vector3::Constant(0.3), Vector3::Constant(0.1))); // 30cm std on x,y,z 0.1 rad on roll,pitch,yaw
- graph.push_back(PriorFactor(Symbol('x', 0), poses[0], poseNoise)); // add directly to graph
+ // 30cm std on x,y,z 0.1 rad on roll,pitch,yaw
+ noiseModel::Diagonal::shared_ptr poseNoise = noiseModel::Diagonal::Sigmas(
+ (Vector(6) << Vector3::Constant(0.3), Vector3::Constant(0.1)));
+ graph.push_back(PriorFactor(0, poses[0], poseNoise));
- // Because the structure-from-motion problem has a scale ambiguity, the problem is still under-constrained
- // Here we add a prior on the second pose x1, so this will fix the scale by indicating the distance between x0 and x1.
- // Because these two are fixed, the rest poses will be alse fixed.
- graph.push_back(PriorFactor(Symbol('x', 1), poses[1], poseNoise)); // add directly to graph
+ // Because the structure-from-motion problem has a scale ambiguity, the problem is
+ // still under-constrained. Here we add a prior on the second pose x1, so this will
+ // fix the scale by indicating the distance between x0 and x1.
+ // Because these two are fixed, the rest of the poses will be also be fixed.
+ graph.push_back(PriorFactor(1, poses[1], poseNoise)); // add directly to graph
graph.print("Factor Graph:\n");
- // Create the data structure to hold the initial estimate to the solution
+ // Create the initial estimate to the solution
// Intentionally initialize the variables off from the ground truth
Values initialEstimate;
+ Pose3 delta(Rot3::rodriguez(-0.1, 0.2, 0.25), Point3(0.05, -0.10, 0.20));
for (size_t i = 0; i < poses.size(); ++i)
- initialEstimate.insert(Symbol('x', i), poses[i].compose(Pose3(Rot3::rodriguez(-0.1, 0.2, 0.25), Point3(0.05, -0.10, 0.20))));
+ initialEstimate.insert(i, poses[i].compose(delta));
initialEstimate.print("Initial Estimates:\n");
// Optimize the graph and print results
Values result = DoglegOptimizer(graph, initialEstimate).optimize();
result.print("Final results:\n");
-
- // Notice: Smart factor represent the 3D point as a factor, not a variable.
+ // A smart factor represent the 3D point inside the factor, not as a variable.
// The 3D position of the landmark is not explicitly calculated by the optimizer.
- // If you do want to output the landmark's 3D position, you should use the internal function point()
- // of the smart factor to get the 3D point.
+ // To obtain the landmark's 3D position, we use the "point" method of the smart factor.
Values landmark_result;
- for (size_t i = 0; i < points.size(); ++i) {
+ for (size_t j = 0; j < points.size(); ++j) {
- // The output of point() is in boost::optional, since sometimes
- // the triangulation opterations inside smart factor will encounter degeneracy.
- // Check the manual of boost::optional for more details for the usages.
+ // The output of point() is in boost::optional, as sometimes
+ // the triangulation operation inside smart factor will encounter degeneracy.
boost::optional point;
- // here we use the saved smart factors to call, pass in all optimized pose to calculate landmark positions
- point = smartfactors_ptr.at(i)->point(result);
-
- // ignore if boost::optional return NULL
- if (point)
- landmark_result.insert(Symbol('l', i), *point);
+ // The graph stores Factor shared_ptrs, so we cast back to a SmartFactor first
+ SmartFactor::shared_ptr smart = boost::dynamic_pointer_cast(graph[j]);
+ if (smart) {
+ point = smart->point(result);
+ if (point) // ignore if boost::optional return NULL
+ landmark_result.insert(j, *point);
+ }
}
landmark_result.print("Landmark results:\n");
-
return 0;
}
/* ************************************************************************* */
diff --git a/examples/SFMExample_bal.cpp b/examples/SFMExample_bal.cpp
index 9d453fe5f..43177dc42 100644
--- a/examples/SFMExample_bal.cpp
+++ b/examples/SFMExample_bal.cpp
@@ -43,8 +43,7 @@ int main (int argc, char* argv[]) {
// Load the SfM data from file
SfM_data mydata;
- const bool success = readBAL(filename, mydata);
- assert(success);
+ assert(readBAL(filename, mydata));
cout << boost::format("read %1% tracks on %2% cameras\n") % mydata.number_tracks() % mydata.number_cameras();
// Create a factor graph
diff --git a/examples/StereoVOExample.cpp b/examples/StereoVOExample.cpp
index 8a04d13a3..19725798c 100644
--- a/examples/StereoVOExample.cpp
+++ b/examples/StereoVOExample.cpp
@@ -24,47 +24,48 @@
*/
#include
-#include
+#include
+#include
+#include
#include
#include
-#include
-#include
-#include
#include
-#include
-#include
-
-
using namespace std;
using namespace gtsam;
int main(int argc, char** argv){
+
//create graph object, add first pose at origin with key '1'
NonlinearFactorGraph graph;
- Pose3 first_pose = Pose3();
- graph.push_back(NonlinearEquality(1, first_pose));
+ Pose3 first_pose;
+ graph.push_back(NonlinearEquality(1, Pose3()));
//create factor noise model with 3 sigmas of value 1
const noiseModel::Isotropic::shared_ptr model = noiseModel::Isotropic::Sigma(3,1);
//create stereo camera calibration object with .2m between cameras
const Cal3_S2Stereo::shared_ptr K(new Cal3_S2Stereo(1000, 1000, 0, 320, 240, 0.2));
+
//create and add stereo factors between first pose (key value 1) and the three landmarks
graph.push_back(GenericStereoFactor(StereoPoint2(520, 480, 440), model, 1, 3, K));
graph.push_back(GenericStereoFactor(StereoPoint2(120, 80, 440), model, 1, 4, K));
graph.push_back(GenericStereoFactor(StereoPoint2(320, 280, 140), model, 1, 5, K));
+
//create and add stereo factors between second pose and the three landmarks
graph.push_back(GenericStereoFactor(StereoPoint2(570, 520, 490), model, 2, 3, K));
graph.push_back(GenericStereoFactor(StereoPoint2(70, 20, 490), model, 2, 4, K));
graph.push_back(GenericStereoFactor(StereoPoint2(320, 270, 115), model, 2, 5, K));
+
//create Values object to contain initial estimates of camera poses and landmark locations
- Values initial_estimate = Values();
+ Values initial_estimate;
+
//create and add iniital estimates
initial_estimate.insert(1, first_pose);
initial_estimate.insert(2, Pose3(Rot3(), Point3(0.1, -0.1, 1.1)));
initial_estimate.insert(3, Point3(1, 1, 5));
initial_estimate.insert(4, Point3(-1, 1, 5));
initial_estimate.insert(5, Point3(0, -0.5, 5));
+
//create Levenberg-Marquardt optimizer for resulting factor graph, optimize
LevenbergMarquardtOptimizer optimizer = LevenbergMarquardtOptimizer(graph, initial_estimate);
Values result = optimizer.optimize();
@@ -72,4 +73,4 @@ int main(int argc, char** argv){
result.print("Final result:\n");
return 0;
-}
\ No newline at end of file
+}
diff --git a/gtsam.h b/gtsam.h
index b7178d534..678e2a3d6 100644
--- a/gtsam.h
+++ b/gtsam.h
@@ -1481,9 +1481,7 @@ class GaussianISAM {
#include
virtual class IterativeOptimizationParameters {
- string getKernel() const ;
string getVerbosity() const;
- void setKernel(string s) ;
void setVerbosity(string s) ;
void print() const;
};
@@ -1516,7 +1514,7 @@ virtual class SubgraphSolverParameters : gtsam::ConjugateGradientParameters {
void print() const;
};
-class SubgraphSolver {
+virtual class SubgraphSolver {
SubgraphSolver(const gtsam::GaussianFactorGraph &A, const gtsam::SubgraphSolverParameters ¶meters, const gtsam::Ordering& ordering);
SubgraphSolver(const gtsam::GaussianFactorGraph &Ab1, const gtsam::GaussianFactorGraph &Ab2, const gtsam::SubgraphSolverParameters ¶meters, const gtsam::Ordering& ordering);
gtsam::VectorValues optimize() const;
@@ -1855,12 +1853,12 @@ virtual class NonlinearOptimizerParams {
void setLinearSolverType(string solver);
void setOrdering(const gtsam::Ordering& ordering);
- void setIterativeParams(const gtsam::SubgraphSolverParameters ¶ms);
+ void setIterativeParams(gtsam::IterativeOptimizationParameters* params);
bool isMultifrontal() const;
bool isSequential() const;
bool isCholmod() const;
- bool isCG() const;
+ bool isIterative() const;
};
bool checkConvergence(double relativeErrorTreshold,
@@ -2199,6 +2197,25 @@ virtual class GeneralSFMFactor2 : gtsam::NoiseModelFactor {
void serialize() const;
};
+#include
+template
+virtual class SmartProjectionPoseFactor : gtsam::NonlinearFactor {
+
+ SmartProjectionPoseFactor(double rankTol, double linThreshold,
+ bool manageDegeneracy, bool enableEPI, const POSE& body_P_sensor);
+
+ SmartProjectionPoseFactor(double rankTol);
+ SmartProjectionPoseFactor();
+
+ void add(const gtsam::Point2& measured_i, size_t poseKey_i, const gtsam::noiseModel::Base* noise_i,
+ const CALIBRATION* K_i);
+
+ // enabling serialization functionality
+ //void serialize() const;
+};
+
+typedef gtsam::SmartProjectionPoseFactor SmartProjectionPose3Factor;
+
#include
template
@@ -2301,7 +2318,8 @@ virtual class ConstantBias : gtsam::Value {
#include
class ImuFactorPreintegratedMeasurements {
// Standard Constructor
- ImuFactorPreintegratedMeasurements(const gtsam::imuBias::ConstantBias& bias, Matrix measuredAccCovariance, Matrix measuredOmegaCovariance, Matrix integrationErrorCovariance);
+ ImuFactorPreintegratedMeasurements(const gtsam::imuBias::ConstantBias& bias, Matrix measuredAccCovariance,Matrix measuredOmegaCovariance, Matrix integrationErrorCovariance, bool use2ndOrderIntegration);
+ ImuFactorPreintegratedMeasurements(const gtsam::imuBias::ConstantBias& bias, Matrix measuredAccCovariance,Matrix measuredOmegaCovariance, Matrix integrationErrorCovariance);
ImuFactorPreintegratedMeasurements(const gtsam::ImuFactorPreintegratedMeasurements& rhs);
// Testable
@@ -2339,6 +2357,15 @@ class CombinedImuFactorPreintegratedMeasurements {
Matrix biasAccCovariance,
Matrix biasOmegaCovariance,
Matrix biasAccOmegaInit);
+ CombinedImuFactorPreintegratedMeasurements(
+ const gtsam::imuBias::ConstantBias& bias,
+ Matrix measuredAccCovariance,
+ Matrix measuredOmegaCovariance,
+ Matrix integrationErrorCovariance,
+ Matrix biasAccCovariance,
+ Matrix biasOmegaCovariance,
+ Matrix biasAccOmegaInit,
+ bool use2ndOrderIntegration);
CombinedImuFactorPreintegratedMeasurements(const gtsam::CombinedImuFactorPreintegratedMeasurements& rhs);
// Testable
@@ -2352,8 +2379,7 @@ class CombinedImuFactorPreintegratedMeasurements {
virtual class CombinedImuFactor : gtsam::NonlinearFactor {
CombinedImuFactor(size_t pose_i, size_t vel_i, size_t pose_j, size_t vel_j, size_t bias_i, size_t bias_j,
- const gtsam::CombinedImuFactorPreintegratedMeasurements& CombinedPreintegratedMeasurements, Vector gravity, Vector omegaCoriolis,
- const gtsam::noiseModel::Base* model);
+ const gtsam::CombinedImuFactorPreintegratedMeasurements& CombinedPreintegratedMeasurements, Vector gravity, Vector omegaCoriolis);
// Standard Interface
gtsam::CombinedImuFactorPreintegratedMeasurements preintegratedMeasurements() const;
diff --git a/gtsam/base/CMakeLists.txt b/gtsam/base/CMakeLists.txt
index 50145846e..66d3ec721 100644
--- a/gtsam/base/CMakeLists.txt
+++ b/gtsam/base/CMakeLists.txt
@@ -7,9 +7,3 @@ install(FILES ${base_headers_tree} DESTINATION include/gtsam/base/treeTraversal)
# Build tests
add_subdirectory(tests)
-
-# Build timing scripts
-if (GTSAM_BUILD_TIMING)
- gtsam_add_subdir_timing(base "gtsam" "gtsam" "${base_excluded_files}")
-endif(GTSAM_BUILD_TIMING)
-
diff --git a/gtsam/discrete/CMakeLists.txt b/gtsam/discrete/CMakeLists.txt
index 3c5602080..d78dff34f 100644
--- a/gtsam/discrete/CMakeLists.txt
+++ b/gtsam/discrete/CMakeLists.txt
@@ -6,9 +6,3 @@ install(FILES ${discrete_headers} DESTINATION include/gtsam/discrete)
# Add all tests
add_subdirectory(tests)
-
-# Build timing scripts
-#if (GTSAM_BUILD_TIMING)
-# gtsam_add_timing(discrete "${discrete_local_libs}")
-#endif(GTSAM_BUILD_TIMING)
-
diff --git a/gtsam/geometry/CMakeLists.txt b/gtsam/geometry/CMakeLists.txt
index f72f965ea..dabdde45c 100644
--- a/gtsam/geometry/CMakeLists.txt
+++ b/gtsam/geometry/CMakeLists.txt
@@ -4,9 +4,3 @@ install(FILES ${geometry_headers} DESTINATION include/gtsam/geometry)
# Build tests
add_subdirectory(tests)
-
-# Build timing scripts
-if (GTSAM_BUILD_TIMING)
- gtsam_add_subdir_timing(geometry "gtsam" "gtsam" "${geometry_excluded_files}")
-endif(GTSAM_BUILD_TIMING)
-
diff --git a/gtsam/geometry/Cal3DS2.cpp b/gtsam/geometry/Cal3DS2.cpp
index e6c97cb1a..c75eff022 100644
--- a/gtsam/geometry/Cal3DS2.cpp
+++ b/gtsam/geometry/Cal3DS2.cpp
@@ -25,7 +25,7 @@ namespace gtsam {
/* ************************************************************************* */
Cal3DS2::Cal3DS2(const Vector &v):
- fx_(v[0]), fy_(v[1]), s_(v[2]), u0_(v[3]), v0_(v[4]), k1_(v[5]), k2_(v[6]), k3_(v[7]), k4_(v[8]){}
+ fx_(v[0]), fy_(v[1]), s_(v[2]), u0_(v[3]), v0_(v[4]), k1_(v[5]), k2_(v[6]), p1_(v[7]), p2_(v[8]){}
/* ************************************************************************* */
Matrix Cal3DS2::K() const {
@@ -34,32 +34,64 @@ Matrix Cal3DS2::K() const {
/* ************************************************************************* */
Vector Cal3DS2::vector() const {
- return (Vector(9) << fx_, fy_, s_, u0_, v0_, k1_, k2_, k3_, k4_);
+ return (Vector(9) << fx_, fy_, s_, u0_, v0_, k1_, k2_, p1_, p2_);
}
/* ************************************************************************* */
-void Cal3DS2::print(const std::string& s) const {
- gtsam::print(K(), s + ".K");
- gtsam::print(Vector(k()), s + ".k");
+void Cal3DS2::print(const std::string& s_) const {
+ gtsam::print(K(), s_ + ".K");
+ gtsam::print(Vector(k()), s_ + ".k");
}
/* ************************************************************************* */
bool Cal3DS2::equals(const Cal3DS2& K, double tol) const {
if (fabs(fx_ - K.fx_) > tol || fabs(fy_ - K.fy_) > tol || fabs(s_ - K.s_) > tol ||
fabs(u0_ - K.u0_) > tol || fabs(v0_ - K.v0_) > tol || fabs(k1_ - K.k1_) > tol ||
- fabs(k2_ - K.k2_) > tol || fabs(k3_ - K.k3_) > tol || fabs(k4_ - K.k4_) > tol)
+ fabs(k2_ - K.k2_) > tol || fabs(p1_ - K.p1_) > tol || fabs(p2_ - K.p2_) > tol)
return false;
return true;
}
/* ************************************************************************* */
-Point2 Cal3DS2::uncalibrate(const Point2& p,
- boost::optional H1,
- boost::optional H2) const {
+static Eigen::Matrix D2dcalibration(double x, double y, double xx,
+ double yy, double xy, double rr, double r4, double pnx, double pny,
+ const Eigen::Matrix& DK) {
+ Eigen::Matrix DR1;
+ DR1 << pnx, 0.0, pny, 1.0, 0.0, 0.0, pny, 0.0, 0.0, 1.0;
+ Eigen::Matrix DR2;
+ DR2 << x * rr, x * r4, 2 * xy, rr + 2 * xx, //
+ y * rr, y * r4, rr + 2 * yy, 2 * xy;
+ Eigen::Matrix D;
+ D << DR1, DK * DR2;
+ return D;
+}
- // parameters
- const double fx = fx_, fy = fy_, s = s_;
- const double k1 = k1_, k2 = k2_, k3 = k3_, k4 = k4_;
+/* ************************************************************************* */
+static Eigen::Matrix D2dintrinsic(double x, double y, double rr,
+ double g, double k1, double k2, double p1, double p2,
+ const Eigen::Matrix& DK) {
+ const double drdx = 2. * x;
+ const double drdy = 2. * y;
+ const double dgdx = k1 * drdx + k2 * 2. * rr * drdx;
+ const double dgdy = k1 * drdy + k2 * 2. * rr * drdy;
+
+ // Dx = 2*p1*xy + p2*(rr+2*xx);
+ // Dy = 2*p2*xy + p1*(rr+2*yy);
+ const double dDxdx = 2. * p1 * y + p2 * (drdx + 4. * x);
+ const double dDxdy = 2. * p1 * x + p2 * drdy;
+ const double dDydx = 2. * p2 * y + p1 * drdx;
+ const double dDydy = 2. * p2 * x + p1 * (drdy + 4. * y);
+
+ Eigen::Matrix DR;
+ DR << g + x * dgdx + dDxdx, x * dgdy + dDxdy, //
+ y * dgdx + dDydx, g + y * dgdy + dDydy;
+
+ return DK * DR;
+}
+
+/* ************************************************************************* */
+Point2 Cal3DS2::uncalibrate(const Point2& p, boost::optional H1,
+ boost::optional H2) const {
// rr = x^2 + y^2;
// g = (1 + k(1)*rr + k(2)*rr^2);
@@ -68,40 +100,29 @@ Point2 Cal3DS2::uncalibrate(const Point2& p,
const double x = p.x(), y = p.y(), xy = x * y, xx = x * x, yy = y * y;
const double rr = xx + yy;
const double r4 = rr * rr;
- const double g = 1. + k1 * rr + k2 * r4;
- const double dx = 2. * k3 * xy + k4 * (rr + 2. * xx);
- const double dy = 2. * k4 * xy + k3 * (rr + 2. * yy);
+ const double g = 1. + k1_ * rr + k2_ * r4; // scaling factor
- const double pnx = g*x + dx;
- const double pny = g*y + dy;
+ // tangential component
+ const double dx = 2. * p1_ * xy + p2_ * (rr + 2. * xx);
+ const double dy = 2. * p2_ * xy + p1_ * (rr + 2. * yy);
- // Inlined derivative for calibration
- if (H1) {
- *H1 = (Matrix(2, 9) << pnx, 0.0, pny, 1.0, 0.0, fx * x * rr + s * y * rr,
- fx * x * r4 + s * y * r4, fx * 2. * xy + s * (rr + 2. * yy),
- fx * (rr + 2. * xx) + s * (2. * xy), 0.0, pny, 0.0, 0.0, 1.0,
- fy * y * rr, fy * y * r4, fy * (rr + 2. * yy), fy * (2. * xy));
- }
- // Inlined derivative for points
- if (H2) {
- const double dr_dx = 2. * x;
- const double dr_dy = 2. * y;
- const double dg_dx = k1 * dr_dx + k2 * 2. * rr * dr_dx;
- const double dg_dy = k1 * dr_dy + k2 * 2. * rr * dr_dy;
+ // Radial and tangential distortion applied
+ const double pnx = g * x + dx;
+ const double pny = g * y + dy;
- const double dDx_dx = 2. * k3 * y + k4 * (dr_dx + 4. * x);
- const double dDx_dy = 2. * k3 * x + k4 * dr_dy;
- const double dDy_dx = 2. * k4 * y + k3 * dr_dx;
- const double dDy_dy = 2. * k4 * x + k3 * (dr_dy + 4. * y);
+ Eigen::Matrix DK;
+ if (H1 || H2) DK << fx_, s_, 0.0, fy_;
- Matrix DK = (Matrix(2, 2) << fx, s_, 0.0, fy);
- Matrix DR = (Matrix(2, 2) << g + x * dg_dx + dDx_dx, x * dg_dy + dDx_dy,
- y * dg_dx + dDy_dx, g + y * dg_dy + dDy_dy);
+ // Derivative for calibration
+ if (H1)
+ *H1 = D2dcalibration(x, y, xx, yy, xy, rr, r4, pnx, pny, DK);
- *H2 = DK * DR;
- }
+ // Derivative for points
+ if (H2)
+ *H2 = D2dintrinsic(x, y, rr, g, k1_, k2_, p1_, p2_, DK);
- return Point2(fx * pnx + s * pny + u0_, fy * pny + v0_);
+ // Regular uncalibrate after distortion
+ return Point2(fx_ * pnx + s_ * pny + u0_, fy_ * pny + v0_);
}
/* ************************************************************************* */
@@ -118,14 +139,14 @@ Point2 Cal3DS2::calibrate(const Point2& pi, const double tol) const {
// iterate until the uncalibrate is close to the actual pixel coordinate
const int maxIterations = 10;
int iteration;
- for ( iteration = 0; iteration < maxIterations; ++iteration ) {
- if ( uncalibrate(pn).distance(pi) <= tol ) break;
- const double x = pn.x(), y = pn.y(), xy = x*y, xx = x*x, yy = y*y;
+ for (iteration = 0; iteration < maxIterations; ++iteration) {
+ if (uncalibrate(pn).distance(pi) <= tol) break;
+ const double x = pn.x(), y = pn.y(), xy = x * y, xx = x * x, yy = y * y;
const double rr = xx + yy;
- const double g = (1+k1_*rr+k2_*rr*rr);
- const double dx = 2*k3_*xy + k4_*(rr+2*xx);
- const double dy = 2*k4_*xy + k3_*(rr+2*yy);
- pn = (invKPi - Point2(dx,dy))/g;
+ const double g = (1 + k1_ * rr + k2_ * rr * rr);
+ const double dx = 2 * p1_ * xy + p2_ * (rr + 2 * xx);
+ const double dy = 2 * p2_ * xy + p1_ * (rr + 2 * yy);
+ pn = (invKPi - Point2(dx, dy)) / g;
}
if ( iteration >= maxIterations )
@@ -136,47 +157,28 @@ Point2 Cal3DS2::calibrate(const Point2& pi, const double tol) const {
/* ************************************************************************* */
Matrix Cal3DS2::D2d_intrinsic(const Point2& p) const {
- //const double fx = fx_, fy = fy_, s = s_;
- const double k1 = k1_, k2 = k2_, k3 = k3_, k4 = k4_;
- //const double x = p.x(), y = p.y(), xx = x*x, yy = y*y, xy = x*y;
- const double x = p.x(), y = p.y(), xx = x*x, yy = y*y;
+ const double x = p.x(), y = p.y(), xx = x * x, yy = y * y;
const double rr = xx + yy;
- const double dr_dx = 2*x;
- const double dr_dy = 2*y;
- const double r4 = rr*rr;
- const double g = 1 + k1*rr + k2*r4;
- const double dg_dx = k1*dr_dx + k2*2*rr*dr_dx;
- const double dg_dy = k1*dr_dy + k2*2*rr*dr_dy;
-
- // Dx = 2*k3*xy + k4*(rr+2*xx);
- // Dy = 2*k4*xy + k3*(rr+2*yy);
- const double dDx_dx = 2*k3*y + k4*(dr_dx + 4*x);
- const double dDx_dy = 2*k3*x + k4*dr_dy;
- const double dDy_dx = 2*k4*y + k3*dr_dx;
- const double dDy_dy = 2*k4*x + k3*(dr_dy + 4*y);
-
- Matrix DK = (Matrix(2, 2) << fx_, s_, 0.0, fy_);
- Matrix DR = (Matrix(2, 2) << g + x*dg_dx + dDx_dx, x*dg_dy + dDx_dy,
- y*dg_dx + dDy_dx, g + y*dg_dy + dDy_dy);
-
- return DK * DR;
+ const double r4 = rr * rr;
+ const double g = (1 + k1_ * rr + k2_ * r4);
+ Eigen::Matrix DK;
+ DK << fx_, s_, 0.0, fy_;
+ return D2dintrinsic(x, y, rr, g, k1_, k2_, p1_, p2_, DK);
}
/* ************************************************************************* */
Matrix Cal3DS2::D2d_calibration(const Point2& p) const {
- const double x = p.x(), y = p.y(), xx = x*x, yy = y*y, xy = x*y;
+ const double x = p.x(), y = p.y(), xx = x * x, yy = y * y, xy = x * y;
const double rr = xx + yy;
- const double r4 = rr*rr;
- const double fx = fx_, fy = fy_, s = s_;
- const double g = (1+k1_*rr+k2_*r4);
- const double dx = 2*k3_*xy + k4_*(rr+2*xx);
- const double dy = 2*k4_*xy + k3_*(rr+2*yy);
- const double pnx = g*x + dx;
- const double pny = g*y + dy;
-
- return (Matrix(2, 9) <<
- pnx, 0.0, pny, 1.0, 0.0, fx*x*rr + s*y*rr, fx*x*r4 + s*y*r4, fx*2*xy + s*(rr+2*yy), fx*(rr+2*xx) + s*(2*xy),
- 0.0, pny, 0.0, 0.0, 1.0, fy*y*rr , fy*y*r4 , fy*(rr+2*yy) , fy*(2*xy) );
+ const double r4 = rr * rr;
+ const double g = (1 + k1_ * rr + k2_ * r4);
+ const double dx = 2 * p1_ * xy + p2_ * (rr + 2 * xx);
+ const double dy = 2 * p2_ * xy + p1_ * (rr + 2 * yy);
+ const double pnx = g * x + dx;
+ const double pny = g * y + dy;
+ Eigen::Matrix DK;
+ DK << fx_, s_, 0.0, fy_;
+ return D2dcalibration(x, y, xx, yy, xy, rr, r4, pnx, pny, DK);
}
/* ************************************************************************* */
diff --git a/gtsam/geometry/Cal3DS2.h b/gtsam/geometry/Cal3DS2.h
index ced05086b..51fe958d6 100644
--- a/gtsam/geometry/Cal3DS2.h
+++ b/gtsam/geometry/Cal3DS2.h
@@ -27,6 +27,15 @@ namespace gtsam {
* @brief Calibration of a camera with radial distortion
* @addtogroup geometry
* \nosubgrouping
+ *
+ * Uses same distortionmodel as OpenCV, with
+ * http://docs.opencv.org/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html
+ * but using only k1,k2,p1, and p2 coefficients.
+ * K = [ fx s u0 ; 0 fy v0 ; 0 0 1 ]
+ * rr = Pn.x^2 + Pn.y^2
+ * \hat{pn} = (1 + k1*rr + k2*rr^2 ) pn + [ 2*k3 pn.x pn.y + k4 (rr + 2 Pn.x^2) ;
+ * k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y ]
+ * pi = K*pn
*/
class GTSAM_EXPORT Cal3DS2 : public DerivedValue {
@@ -34,28 +43,22 @@ protected:
double fx_, fy_, s_, u0_, v0_ ; // focal length, skew and principal point
double k1_, k2_ ; // radial 2nd-order and 4th-order
- double k3_, k4_ ; // tangential distortion
-
- // K = [ fx s u0 ; 0 fy v0 ; 0 0 1 ]
- // rr = Pn.x^2 + Pn.y^2
- // \hat{pn} = (1 + k1*rr + k2*rr^2 ) pn + [ 2*k3 pn.x pn.y + k4 (rr + 2 Pn.x^2) ;
- // k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y ]
- // pi = K*pn
+ double p1_, p2_ ; // tangential distortion
public:
Matrix K() const ;
- Eigen::Vector4d k() const { return Eigen::Vector4d(k1_, k2_, k3_, k4_); }
+ Eigen::Vector4d k() const { return Eigen::Vector4d(k1_, k2_, p1_, p2_); }
Vector vector() const ;
/// @name Standard Constructors
/// @{
/// Default Constructor with only unit focal length
- Cal3DS2() : fx_(1), fy_(1), s_(0), u0_(0), v0_(0), k1_(0), k2_(0), k3_(0), k4_(0) {}
+ Cal3DS2() : fx_(1), fy_(1), s_(0), u0_(0), v0_(0), k1_(0), k2_(0), p1_(0), p2_(0) {}
Cal3DS2(double fx, double fy, double s, double u0, double v0,
- double k1, double k2, double k3 = 0.0, double k4 = 0.0) :
- fx_(fx), fy_(fy), s_(s), u0_(u0), v0_(v0), k1_(k1), k2_(k2), k3_(k3), k4_(k4) {}
+ double k1, double k2, double p1 = 0.0, double p2 = 0.0) :
+ fx_(fx), fy_(fy), s_(s), u0_(u0), v0_(v0), k1_(k1), k2_(k2), p1_(p1), p2_(p2) {}
/// @}
/// @name Advanced Constructors
@@ -92,18 +95,30 @@ public:
/// image center in y
inline double py() const { return v0_;}
+ /// First distortion coefficient
+ inline double k1() const { return k1_;}
+
+ /// Second distortion coefficient
+ inline double k2() const { return k2_;}
+
+ /// First tangential distortion coefficient
+ inline double p1() const { return p1_;}
+
+ /// Second tangential distortion coefficient
+ inline double p2() const { return p2_;}
+
/**
- * convert intrinsic coordinates xy to image coordinates uv
+ * convert intrinsic coordinates xy to (distorted) image coordinates uv
* @param p point in intrinsic coordinates
* @param Dcal optional 2*9 Jacobian wrpt Cal3DS2 parameters
* @param Dp optional 2*2 Jacobian wrpt intrinsic coordinates
- * @return point in image coordinates
+ * @return point in (distorted) image coordinates
*/
Point2 uncalibrate(const Point2& p,
boost::optional Dcal = boost::none,
boost::optional Dp = boost::none) const ;
- /// Conver a pixel coordinate to ideal coordinate
+ /// Convert (distorted) image coordinates uv to intrinsic coordinates xy
Point2 calibrate(const Point2& p, const double tol=1e-5) const;
/// Derivative of uncalibrate wrpt intrinsic coordinates
@@ -148,8 +163,8 @@ private:
ar & BOOST_SERIALIZATION_NVP(v0_);
ar & BOOST_SERIALIZATION_NVP(k1_);
ar & BOOST_SERIALIZATION_NVP(k2_);
- ar & BOOST_SERIALIZATION_NVP(k3_);
- ar & BOOST_SERIALIZATION_NVP(k4_);
+ ar & BOOST_SERIALIZATION_NVP(p1_);
+ ar & BOOST_SERIALIZATION_NVP(p2_);
}
diff --git a/gtsam/geometry/Cal3Unified.cpp b/gtsam/geometry/Cal3Unified.cpp
index 808cb84a4..e7b408982 100644
--- a/gtsam/geometry/Cal3Unified.cpp
+++ b/gtsam/geometry/Cal3Unified.cpp
@@ -43,7 +43,7 @@ void Cal3Unified::print(const std::string& s) const {
bool Cal3Unified::equals(const Cal3Unified& K, double tol) const {
if (fabs(fx_ - K.fx_) > tol || fabs(fy_ - K.fy_) > tol || fabs(s_ - K.s_) > tol ||
fabs(u0_ - K.u0_) > tol || fabs(v0_ - K.v0_) > tol || fabs(k1_ - K.k1_) > tol ||
- fabs(k2_ - K.k2_) > tol || fabs(k3_ - K.k3_) > tol || fabs(k4_ - K.k4_) > tol ||
+ fabs(k2_ - K.k2_) > tol || fabs(p1_ - K.p1_) > tol || fabs(p2_ - K.p2_) > tol ||
fabs(xi_ - K.xi_) > tol)
return false;
return true;
diff --git a/gtsam/geometry/Cal3Unified.h b/gtsam/geometry/Cal3Unified.h
index a1d47332a..58e024c27 100644
--- a/gtsam/geometry/Cal3Unified.h
+++ b/gtsam/geometry/Cal3Unified.h
@@ -31,6 +31,14 @@ namespace gtsam {
* @brief Calibration of a omni-directional camera with mirror + lens radial distortion
* @addtogroup geometry
* \nosubgrouping
+ *
+ * Similar to Cal3DS2, does distortion but has additional mirror parameter xi
+ * K = [ fx s u0 ; 0 fy v0 ; 0 0 1 ]
+ * Pn = [ P.x / (1 + xi * \sqrt{P.x^2 + P.y^2 + 1}), P.y / (1 + xi * \sqrt{P.x^2 + P.y^2 + 1})]
+ * rr = Pn.x^2 + Pn.y^2
+ * \hat{pn} = (1 + k1*rr + k2*rr^2 ) pn + [ 2*k3 pn.x pn.y + k4 (rr + 2 Pn.x^2) ;
+ * k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y ]
+ * pi = K*pn
*/
class GTSAM_EXPORT Cal3Unified : public Cal3DS2 {
@@ -41,13 +49,6 @@ private:
double xi_; // mirror parameter
- // K = [ fx s u0 ; 0 fy v0 ; 0 0 1 ]
- // Pn = [ P.x / (1 + xi * \sqrt{P.x^2 + P.y^2 + 1}), P.y / (1 + xi * \sqrt{P.x^2 + P.y^2 + 1})]
- // rr = Pn.x^2 + Pn.y^2
- // \hat{pn} = (1 + k1*rr + k2*rr^2 ) pn + [ 2*k3 pn.x pn.y + k4 (rr + 2 Pn.x^2) ;
- // k3 (rr + 2 Pn.y^2) + 2*k4 pn.x pn.y ]
- // pi = K*pn
-
public:
//Matrix K() const ;
//Eigen::Vector4d k() const { return Base::k(); }
@@ -60,8 +61,8 @@ public:
Cal3Unified() : Base(), xi_(0) {}
Cal3Unified(double fx, double fy, double s, double u0, double v0,
- double k1, double k2, double k3 = 0.0, double k4 = 0.0, double xi = 0.0) :
- Base(fx, fy, s, u0, v0, k1, k2, k3, k4), xi_(xi) {}
+ double k1, double k2, double p1 = 0.0, double p2 = 0.0, double xi = 0.0) :
+ Base(fx, fy, s, u0, v0, k1, k2, p1, p2), xi_(xi) {}
/// @}
/// @name Advanced Constructors
diff --git a/gtsam/geometry/EssentialMatrix.h b/gtsam/geometry/EssentialMatrix.h
index 3c9568c3d..32b966261 100644
--- a/gtsam/geometry/EssentialMatrix.h
+++ b/gtsam/geometry/EssentialMatrix.h
@@ -112,6 +112,16 @@ public:
return E_;
}
+ /// Return epipole in image_a , as Unit3 to allow for infinity
+ inline const Unit3& epipole_a() const {
+ return aTb_; // == direction()
+ }
+
+ /// Return epipole in image_b, as Unit3 to allow for infinity
+ inline Unit3 epipole_b() const {
+ return aRb_.unrotate(aTb_); // == rotation.unrotate(direction())
+ }
+
/**
* @brief takes point in world coordinates and transforms it to pose with |t|==1
* @param p point in world coordinates
diff --git a/gtsam/geometry/PinholeCamera.h b/gtsam/geometry/PinholeCamera.h
index 19eb87b5f..709b95181 100644
--- a/gtsam/geometry/PinholeCamera.h
+++ b/gtsam/geometry/PinholeCamera.h
@@ -275,11 +275,12 @@ public:
if (P.z() <= 0)
throw CheiralityException();
#endif
+ double d = 1.0 / P.z();
+ const double u = P.x() * d, v = P.y() * d;
if (Dpoint) {
- double d = 1.0 / P.z(), d2 = d * d;
- *Dpoint = (Matrix(2, 3) << d, 0.0, -P.x() * d2, 0.0, d, -P.y() * d2);
+ *Dpoint = (Matrix(2, 3) << d, 0.0, -u * d, 0.0, d, -v * d);
}
- return Point2(P.x() / P.z(), P.y() / P.z());
+ return Point2(u, v);
}
/// Project a point into the image and check depth
diff --git a/gtsam/geometry/Unit3.cpp b/gtsam/geometry/Unit3.cpp
index db4c8da83..ad6a41bde 100644
--- a/gtsam/geometry/Unit3.cpp
+++ b/gtsam/geometry/Unit3.cpp
@@ -58,7 +58,7 @@ Unit3 Unit3::Random(boost::mt19937 & rng) {
}
/* ************************************************************************* */
-Matrix Unit3::basis() const {
+const Matrix& Unit3::basis() const {
// Return cached version if exists
if (B_.rows() == 3)
diff --git a/gtsam/geometry/Unit3.h b/gtsam/geometry/Unit3.h
index 39e2c9799..8d2c024c0 100644
--- a/gtsam/geometry/Unit3.h
+++ b/gtsam/geometry/Unit3.h
@@ -84,7 +84,7 @@ public:
* It is a 3*2 matrix [b1 b2] composed of two orthogonal directions
* tangent to the sphere at the current direction.
*/
- Matrix basis() const;
+ const Matrix& basis() const;
/// Return skew-symmetric associated with 3D point on unit sphere
Matrix skew() const;
diff --git a/gtsam/geometry/tests/testCal3DS2.cpp b/gtsam/geometry/tests/testCal3DS2.cpp
index fc457767e..eddd3a674 100644
--- a/gtsam/geometry/tests/testCal3DS2.cpp
+++ b/gtsam/geometry/tests/testCal3DS2.cpp
@@ -60,6 +60,8 @@ TEST( Cal3DS2, Duncalibrate1)
K.uncalibrate(p, computed, boost::none);
Matrix numerical = numericalDerivative21(uncalibrate_, K, p, 1e-7);
CHECK(assert_equal(numerical,computed,1e-5));
+ Matrix separate = K.D2d_calibration(p);
+ CHECK(assert_equal(numerical,separate,1e-5));
}
/* ************************************************************************* */
@@ -68,6 +70,8 @@ TEST( Cal3DS2, Duncalibrate2)
Matrix computed; K.uncalibrate(p, boost::none, computed);
Matrix numerical = numericalDerivative22(uncalibrate_, K, p, 1e-7);
CHECK(assert_equal(numerical,computed,1e-5));
+ Matrix separate = K.D2d_intrinsic(p);
+ CHECK(assert_equal(numerical,separate,1e-5));
}
/* ************************************************************************* */
diff --git a/gtsam/geometry/tests/testEssentialMatrix.cpp b/gtsam/geometry/tests/testEssentialMatrix.cpp
index 9ad30bc41..4d34bbe3d 100644
--- a/gtsam/geometry/tests/testEssentialMatrix.cpp
+++ b/gtsam/geometry/tests/testEssentialMatrix.cpp
@@ -144,9 +144,9 @@ TEST (EssentialMatrix, FromPose3_b) {
Matrix actualH;
Rot3 c1Rc2 = Rot3::ypr(0.1, -0.2, 0.3);
Point3 c1Tc2(0.4, 0.5, 0.6);
- EssentialMatrix trueE(c1Rc2, Unit3(c1Tc2));
+ EssentialMatrix E(c1Rc2, Unit3(c1Tc2));
Pose3 pose(c1Rc2, c1Tc2); // Pose between two cameras
- EXPECT(assert_equal(trueE, EssentialMatrix::FromPose3(pose, actualH), 1e-8));
+ EXPECT(assert_equal(E, EssentialMatrix::FromPose3(pose, actualH), 1e-8));
Matrix expectedH = numericalDerivative11(
boost::bind(EssentialMatrix::FromPose3, _1, boost::none), pose);
EXPECT(assert_equal(expectedH, actualH, 1e-5));
@@ -161,6 +161,35 @@ TEST (EssentialMatrix, streaming) {
EXPECT(assert_equal(expected, actual));
}
+//*************************************************************************
+TEST (EssentialMatrix, epipoles) {
+ // Create an E
+ Rot3 c1Rc2 = Rot3::ypr(0.1, -0.2, 0.3);
+ Point3 c1Tc2(0.4, 0.5, 0.6);
+ EssentialMatrix E(c1Rc2, Unit3(c1Tc2));
+
+ // Calculate expected values through SVD
+ Matrix U, V;
+ Vector S;
+ gtsam::svd(E.matrix(), U, S, V);
+
+ // check rank 2 constraint
+ CHECK(fabs(S(2))<1e-10);
+
+ // check epipoles not at infinity
+ CHECK(fabs(U(2,2))>1e-10 && fabs(V(2,2))>1e-10);
+
+ // Check epipoles
+
+ // Epipole in image 1 is just E.direction()
+ Unit3 e1(U(0, 2), U(1, 2), U(2, 2));
+ EXPECT(assert_equal(e1, E.epipole_a()));
+
+ // Epipole in image 2 is E.rotation().unrotate(E.direction())
+ Unit3 e2(V(0, 2), V(1, 2), V(2, 2));
+ EXPECT(assert_equal(e2, E.epipole_b()));
+}
+
/* ************************************************************************* */
int main() {
TestResult tr;
diff --git a/gtsam/geometry/tests/testUnit3.cpp b/gtsam/geometry/tests/testUnit3.cpp
index eb45ea60f..f8b3e7fa4 100644
--- a/gtsam/geometry/tests/testUnit3.cpp
+++ b/gtsam/geometry/tests/testUnit3.cpp
@@ -26,7 +26,7 @@
#include
#include
#include
-#include
+//#include
#include
#include
@@ -234,7 +234,7 @@ TEST(Unit3, localCoordinates_retract) {
for (size_t i = 0; i < numIterations; i++) {
// Sleep for the random number generator (TODO?: Better create all of them first).
- boost::this_thread::sleep(boost::posix_time::milliseconds(0));
+// boost::this_thread::sleep(boost::posix_time::milliseconds(0));
// Create the two Unit3s.
// NOTE: You can not create two totally random Unit3's because you cannot always compute
@@ -264,10 +264,10 @@ TEST(Unit3, localCoordinates_retract_expmap) {
for (size_t i = 0; i < numIterations; i++) {
// Sleep for the random number generator (TODO?: Better create all of them first).
- boost::this_thread::sleep(boost::posix_time::milliseconds(0));
+// boost::this_thread::sleep(boost::posix_time::milliseconds(0));
// Create the two Unit3s.
- // Unlike the above case, we can use any two sphers.
+ // Unlike the above case, we can use any two Unit3's.
Unit3 s1(Point3(randomVector(minSphereLimit, maxSphereLimit)));
// Unit3 s2 (Point3(randomVector(minSphereLimit, maxSphereLimit)));
Vector v12 = randomVector(minXiLimit, maxXiLimit);
diff --git a/gtsam/inference/CMakeLists.txt b/gtsam/inference/CMakeLists.txt
index d5c37d976..c3df3a989 100644
--- a/gtsam/inference/CMakeLists.txt
+++ b/gtsam/inference/CMakeLists.txt
@@ -4,9 +4,3 @@ install(FILES ${inference_headers} DESTINATION include/gtsam/inference)
# Build tests
add_subdirectory(tests)
-
-# Build timing scripts
-if (GTSAM_BUILD_TIMING)
- gtsam_add_subdir_timing(inference "gtsam" "gtsam" "${inference_excluded_files}")
-endif(GTSAM_BUILD_TIMING)
-
diff --git a/gtsam/inference/Ordering.h b/gtsam/inference/Ordering.h
index ea2ad7eda..7b1a2bb2e 100644
--- a/gtsam/inference/Ordering.h
+++ b/gtsam/inference/Ordering.h
@@ -17,9 +17,11 @@
#pragma once
+#include
#include
#include
+#include
#include
#include
#include
@@ -135,6 +137,15 @@ namespace gtsam {
static GTSAM_EXPORT Ordering COLAMDConstrained(const VariableIndex& variableIndex,
const FastMap& groups);
+ /// Return a natural Ordering. Typically used by iterative solvers
+ template
+ static Ordering Natural(const FactorGraph &fg) {
+ FastSet src = fg.keys();
+ std::vector keys(src.begin(), src.end());
+ std::stable_sort(keys.begin(), keys.end());
+ return Ordering(keys);
+ }
+
/// @}
/// @name Testable @{
diff --git a/gtsam/linear/CMakeLists.txt b/gtsam/linear/CMakeLists.txt
index 29eb60b93..084c27057 100644
--- a/gtsam/linear/CMakeLists.txt
+++ b/gtsam/linear/CMakeLists.txt
@@ -4,8 +4,3 @@ install(FILES ${linear_headers} DESTINATION include/gtsam/linear)
# Build tests
add_subdirectory(tests)
-
-# Build timing scripts
-if (GTSAM_BUILD_TIMING)
- gtsam_add_subdir_timing(linear "gtsam" "gtsam" "${linear_excluded_files}")
-endif(GTSAM_BUILD_TIMING)
diff --git a/gtsam/linear/ConjugateGradientSolver.cpp b/gtsam/linear/ConjugateGradientSolver.cpp
new file mode 100644
index 000000000..0505f6c01
--- /dev/null
+++ b/gtsam/linear/ConjugateGradientSolver.cpp
@@ -0,0 +1,49 @@
+/*
+ * ConjugateGradientSolver.cpp
+ *
+ * Created on: Jun 4, 2014
+ * Author: ydjian
+ */
+
+#include
+#include
+#include
+
+using namespace std;
+
+namespace gtsam {
+
+/*****************************************************************************/
+void ConjugateGradientParameters::print(ostream &os) const {
+ Base::print(os);
+ cout << "ConjugateGradientParameters" << endl
+ << "minIter: " << minIterations_ << endl
+ << "maxIter: " << maxIterations_ << endl
+ << "resetIter: " << reset_ << endl
+ << "eps_rel: " << epsilon_rel_ << endl
+ << "eps_abs: " << epsilon_abs_ << endl;
+}
+
+/*****************************************************************************/
+std::string ConjugateGradientParameters::blasTranslator(const BLASKernel value) {
+ std::string s;
+ switch (value) {
+ case ConjugateGradientParameters::GTSAM: s = "GTSAM" ; break;
+ default: s = "UNDEFINED" ; break;
+ }
+ return s;
+}
+
+/*****************************************************************************/
+ConjugateGradientParameters::BLASKernel ConjugateGradientParameters::blasTranslator(const std::string &src) {
+ std::string s = src; boost::algorithm::to_upper(s);
+ if (s == "GTSAM") return ConjugateGradientParameters::GTSAM;
+
+ /* default is SBM */
+ return ConjugateGradientParameters::GTSAM;
+}
+
+
+}
+
+
diff --git a/gtsam/linear/ConjugateGradientSolver.h b/gtsam/linear/ConjugateGradientSolver.h
index d1b3b2c7e..6e8509309 100644
--- a/gtsam/linear/ConjugateGradientSolver.h
+++ b/gtsam/linear/ConjugateGradientSolver.h
@@ -12,14 +12,15 @@
#pragma once
#include
+#include
namespace gtsam {
/**
* parameters for the conjugate gradient method
*/
+class GTSAM_EXPORT ConjugateGradientParameters : public IterativeOptimizationParameters {
-class ConjugateGradientParameters : public IterativeOptimizationParameters {
public:
typedef IterativeOptimizationParameters Base;
typedef boost::shared_ptr shared_ptr;
@@ -30,14 +31,23 @@ public:
double epsilon_rel_; ///< threshold for relative error decrease
double epsilon_abs_; ///< threshold for absolute error decrease
- ConjugateGradientParameters()
- : minIterations_(1), maxIterations_(500), reset_(501), epsilon_rel_(1e-3), epsilon_abs_(1e-3){}
+ /* Matrix Operation Kernel */
+ enum BLASKernel {
+ GTSAM = 0, ///< Jacobian Factor Graph of GTSAM
+ } blas_kernel_ ;
- ConjugateGradientParameters(size_t minIterations, size_t maxIterations, size_t reset, double epsilon_rel, double epsilon_abs)
- : minIterations_(minIterations), maxIterations_(maxIterations), reset_(reset), epsilon_rel_(epsilon_rel), epsilon_abs_(epsilon_abs){}
+ ConjugateGradientParameters()
+ : minIterations_(1), maxIterations_(500), reset_(501), epsilon_rel_(1e-3),
+ epsilon_abs_(1e-3), blas_kernel_(GTSAM) {}
+
+ ConjugateGradientParameters(size_t minIterations, size_t maxIterations, size_t reset,
+ double epsilon_rel, double epsilon_abs, BLASKernel blas)
+ : minIterations_(minIterations), maxIterations_(maxIterations), reset_(reset),
+ epsilon_rel_(epsilon_rel), epsilon_abs_(epsilon_abs), blas_kernel_(blas) {}
ConjugateGradientParameters(const ConjugateGradientParameters &p)
- : Base(p), minIterations_(p.minIterations_), maxIterations_(p.maxIterations_), reset_(p.reset_), epsilon_rel_(p.epsilon_rel_), epsilon_abs_(p.epsilon_abs_) {}
+ : Base(p), minIterations_(p.minIterations_), maxIterations_(p.maxIterations_), reset_(p.reset_),
+ epsilon_rel_(p.epsilon_rel_), epsilon_abs_(p.epsilon_abs_), blas_kernel_(GTSAM) {}
/* general interface */
inline size_t minIterations() const { return minIterations_; }
@@ -61,15 +71,81 @@ public:
inline void setEpsilon_rel(double value) { epsilon_rel_ = value; }
inline void setEpsilon_abs(double value) { epsilon_abs_ = value; }
- virtual void print() const {
- Base::print();
- std::cout << "ConjugateGradientParameters" << std::endl
- << "minIter: " << minIterations_ << std::endl
- << "maxIter: " << maxIterations_ << std::endl
- << "resetIter: " << reset_ << std::endl
- << "eps_rel: " << epsilon_rel_ << std::endl
- << "eps_abs: " << epsilon_abs_ << std::endl;
- }
+
+ void print() const { Base::print(); }
+ virtual void print(std::ostream &os) const;
+
+ static std::string blasTranslator(const BLASKernel k) ;
+ static BLASKernel blasTranslator(const std::string &s) ;
};
+/*********************************************************************************************/
+/*
+ * A template of linear preconditioned conjugate gradient method.
+ * System class should support residual(v, g), multiply(v,Av), leftPrecondition(v, S^{-t}v,
+ * rightPrecondition(v, S^{-1}v), scal(alpha,v), dot(v,v), axpy(alpha,x,y)
+ * Note that the residual is in the preconditioned domain. Refer to Section 9.2 of Saad's book.
+ */
+template
+V preconditionedConjugateGradient(const S &system, const V &initial, const ConjugateGradientParameters ¶meters) {
+
+ V estimate, residual, direction, q1, q2;
+ estimate = residual = direction = q1 = q2 = initial;
+
+ system.residual(estimate, q1); /* q1 = b-Ax */
+ system.leftPrecondition(q1, residual); /* r = S^{-T} (b-Ax) */
+ system.rightPrecondition(residual, direction);/* d = S^{-1} r */
+
+ double currentGamma = system.dot(residual, residual), prevGamma, alpha, beta;
+
+ const size_t iMaxIterations = parameters.maxIterations(),
+ iMinIterations = parameters.minIterations(),
+ iReset = parameters.reset() ;
+ const double threshold = std::max(parameters.epsilon_abs(),
+ parameters.epsilon() * parameters.epsilon() * currentGamma);
+
+ if (parameters.verbosity() >= ConjugateGradientParameters::COMPLEXITY )
+ std::cout << "[PCG] epsilon = " << parameters.epsilon()
+ << ", max = " << parameters.maxIterations()
+ << ", reset = " << parameters.reset()
+ << ", ||r0||^2 = " << currentGamma
+ << ", threshold = " << threshold << std::endl;
+
+ size_t k;
+ for ( k = 1 ; k <= iMaxIterations && (currentGamma > threshold || k <= iMinIterations) ; k++ ) {
+
+ if ( k % iReset == 0 ) {
+ system.residual(estimate, q1); /* q1 = b-Ax */
+ system.leftPrecondition(q1, residual); /* r = S^{-T} (b-Ax) */
+ system.rightPrecondition(residual, direction); /* d = S^{-1} r */
+ currentGamma = system.dot(residual, residual);
+ }
+ system.multiply(direction, q1); /* q1 = A d */
+ alpha = currentGamma / system.dot(direction, q1); /* alpha = gamma / (d' A d) */
+ system.axpy(alpha, direction, estimate); /* estimate += alpha * direction */
+ system.leftPrecondition(q1, q2); /* q2 = S^{-T} * q1 */
+ system.axpy(-alpha, q2, residual); /* residual -= alpha * q2 */
+ prevGamma = currentGamma;
+ currentGamma = system.dot(residual, residual); /* gamma = |residual|^2 */
+ beta = currentGamma / prevGamma;
+ system.rightPrecondition(residual, q1); /* q1 = S^{-1} residual */
+ system.scal(beta, direction);
+ system.axpy(1.0, q1, direction); /* direction = q1 + beta * direction */
+
+ if (parameters.verbosity() >= ConjugateGradientParameters::ERROR )
+ std::cout << "[PCG] k = " << k
+ << ", alpha = " << alpha
+ << ", beta = " << beta
+ << ", ||r||^2 = " << currentGamma
+ << std::endl;
+ }
+ if (parameters.verbosity() >= ConjugateGradientParameters::COMPLEXITY )
+ std::cout << "[PCG] iterations = " << k
+ << ", ||r||^2 = " << currentGamma
+ << std::endl;
+
+ return estimate;
+}
+
+
}
diff --git a/gtsam/linear/GaussianFactorGraph.cpp b/gtsam/linear/GaussianFactorGraph.cpp
index d6836783b..438a06783 100644
--- a/gtsam/linear/GaussianFactorGraph.cpp
+++ b/gtsam/linear/GaussianFactorGraph.cpp
@@ -54,6 +54,20 @@ namespace gtsam {
return keys;
}
+ /* ************************************************************************* */
+ std::map GaussianFactorGraph::getKeyDimMap() const {
+ map spec;
+ BOOST_FOREACH ( const GaussianFactor::shared_ptr &gf, *this ) {
+ for ( GaussianFactor::const_iterator it = gf->begin() ; it != gf->end() ; it++ ) {
+ map::iterator it2 = spec.find(*it);
+ if ( it2 == spec.end() ) {
+ spec.insert(make_pair(*it, gf->getDim(it)));
+ }
+ }
+ }
+ return spec;
+ }
+
/* ************************************************************************* */
vector GaussianFactorGraph::getkeydim() const {
// First find dimensions of each variable
diff --git a/gtsam/linear/GaussianFactorGraph.h b/gtsam/linear/GaussianFactorGraph.h
index 692310f26..910b25d1e 100644
--- a/gtsam/linear/GaussianFactorGraph.h
+++ b/gtsam/linear/GaussianFactorGraph.h
@@ -138,6 +138,9 @@ namespace gtsam {
typedef FastSet Keys;
Keys keys() const;
+ /* return a map of (Key, dimension) */
+ std::map getKeyDimMap() const;
+
std::vector getkeydim() const;
/** unnormalized error */
diff --git a/gtsam/linear/IterativeSolver.cpp b/gtsam/linear/IterativeSolver.cpp
index 2bba3e12b..ab3ccde13 100644
--- a/gtsam/linear/IterativeSolver.cpp
+++ b/gtsam/linear/IterativeSolver.cpp
@@ -6,25 +6,42 @@
*/
#include
+#include
+#include
#include
+#include
#include
+using namespace std;
+
namespace gtsam {
-std::string IterativeOptimizationParameters::getKernel() const { return kernelTranslator(kernel_); }
-std::string IterativeOptimizationParameters::getVerbosity() const { return verbosityTranslator(verbosity_); }
-void IterativeOptimizationParameters::setKernel(const std::string &src) { kernel_ = kernelTranslator(src); }
-void IterativeOptimizationParameters::setVerbosity(const std::string &src) { verbosity_ = verbosityTranslator(src); }
+/*****************************************************************************/
+string IterativeOptimizationParameters::getVerbosity() const { return verbosityTranslator(verbosity_); }
-IterativeOptimizationParameters::Kernel IterativeOptimizationParameters::kernelTranslator(const std::string &src) {
- std::string s = src; boost::algorithm::to_upper(s);
- if (s == "CG") return IterativeOptimizationParameters::CG;
- /* default is cg */
- else return IterativeOptimizationParameters::CG;
+/*****************************************************************************/
+void IterativeOptimizationParameters::setVerbosity(const string &src) { verbosity_ = verbosityTranslator(src); }
+
+/*****************************************************************************/
+void IterativeOptimizationParameters::print() const {
+ print(cout);
}
-IterativeOptimizationParameters::Verbosity IterativeOptimizationParameters::verbosityTranslator(const std::string &src) {
- std::string s = src; boost::algorithm::to_upper(s);
+/*****************************************************************************/
+void IterativeOptimizationParameters::print(ostream &os) const {
+ os << "IterativeOptimizationParameters:" << endl
+ << "verbosity: " << verbosityTranslator(verbosity_) << endl;
+}
+
+/*****************************************************************************/
+ ostream& operator<<(ostream &os, const IterativeOptimizationParameters &p) {
+ p.print(os);
+ return os;
+}
+
+/*****************************************************************************/
+IterativeOptimizationParameters::Verbosity IterativeOptimizationParameters::verbosityTranslator(const string &src) {
+ string s = src; boost::algorithm::to_upper(s);
if (s == "SILENT") return IterativeOptimizationParameters::SILENT;
else if (s == "COMPLEXITY") return IterativeOptimizationParameters::COMPLEXITY;
else if (s == "ERROR") return IterativeOptimizationParameters::ERROR;
@@ -32,18 +49,85 @@ IterativeOptimizationParameters::Verbosity IterativeOptimizationParameters::verb
else return IterativeOptimizationParameters::SILENT;
}
-std::string IterativeOptimizationParameters::kernelTranslator(IterativeOptimizationParameters::Kernel k) {
- if ( k == CG ) return "CG";
- else return "UNKNOWN";
-}
-
-std::string IterativeOptimizationParameters::verbosityTranslator(IterativeOptimizationParameters::Verbosity verbosity) {
+/*****************************************************************************/
+ string IterativeOptimizationParameters::verbosityTranslator(IterativeOptimizationParameters::Verbosity verbosity) {
if (verbosity == SILENT) return "SILENT";
else if (verbosity == COMPLEXITY) return "COMPLEXITY";
else if (verbosity == ERROR) return "ERROR";
else return "UNKNOWN";
}
+/*****************************************************************************/
+VectorValues IterativeSolver::optimize(
+ const GaussianFactorGraph &gfg,
+ boost::optional keyInfo,
+ boost::optional&> lambda)
+{
+ return optimize(
+ gfg,
+ keyInfo ? *keyInfo : KeyInfo(gfg),
+ lambda ? *lambda : std::map());
+}
+
+/*****************************************************************************/
+VectorValues IterativeSolver::optimize (
+ const GaussianFactorGraph &gfg,
+ const KeyInfo &keyInfo,
+ const std::map &lambda)
+{
+ return optimize(gfg, keyInfo, lambda, keyInfo.x0());
+}
+
+/****************************************************************************/
+KeyInfo::KeyInfo(const GaussianFactorGraph &fg, const Ordering &ordering)
+ : ordering_(ordering) {
+ initialize(fg);
+}
+
+/****************************************************************************/
+KeyInfo::KeyInfo(const GaussianFactorGraph &fg)
+ : ordering_(Ordering::Natural(fg)) {
+ initialize(fg);
+}
+
+/****************************************************************************/
+void KeyInfo::initialize(const GaussianFactorGraph &fg){
+ const map colspec = fg.getKeyDimMap();
+ const size_t n = ordering_.size();
+ size_t start = 0;
+
+ for ( size_t i = 0 ; i < n ; ++i ) {
+ const size_t key = ordering_[i];
+ const size_t dim = colspec.find(key)->second;
+ insert(make_pair(key, KeyInfoEntry(i, dim, start)));
+ start += dim ;
+ }
+ numCols_ = start;
+}
+
+/****************************************************************************/
+vector KeyInfo::colSpec() const {
+ std::vector result(size(), 0);
+ BOOST_FOREACH ( const gtsam::KeyInfo::value_type &item, *this ) {
+ result[item.second.index()] = item.second.dim();
+ }
+ return result;
+}
+
+/****************************************************************************/
+VectorValues KeyInfo::x0() const {
+ VectorValues result;
+ BOOST_FOREACH ( const gtsam::KeyInfo::value_type &item, *this ) {
+ result.insert(item.first, Vector::Zero(item.second.dim()));
+ }
+ return result;
+}
+
+/****************************************************************************/
+Vector KeyInfo::x0vector() const {
+ return Vector::Zero(numCols_);
+}
+
}
diff --git a/gtsam/linear/IterativeSolver.h b/gtsam/linear/IterativeSolver.h
index 988293a4f..3c397c9e9 100644
--- a/gtsam/linear/IterativeSolver.h
+++ b/gtsam/linear/IterativeSolver.h
@@ -11,17 +11,27 @@
#pragma once
+#include
#include
-
+#include
+#include
+#include
+#include
+#include
+#include