diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0238ce8 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +.vscode/* +imgui.ini +rviz/* diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..a4b66c0 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,166 @@ +# SPDX-License-Identifier: BSD-2-Clause +cmake_minimum_required(VERSION 2.8.3) +project(hdl_graph_slam) + +# Can we use C++17 in indigo? +add_definitions(-msse -msse2 -msse3 -msse4 -msse4.1 -msse4.2) +set(CMAKE_CXX_FLAGS "-msse -msse2 -msse3 -msse4 -msse4.1 -msse4.2") +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_LIST_DIR}/cmake") + +# pcl 1.7 causes a segfault when it is built with debug mode +set(CMAKE_BUILD_TYPE "RELEASE") + +find_package(catkin REQUIRED COMPONENTS + roscpp + rospy + pcl_ros + geodesy + nmea_msgs + sensor_msgs + geometry_msgs + message_generation + interactive_markers + ndt_omp + fast_gicp +) +catkin_python_setup() + +find_package(PCL REQUIRED) +include_directories(${PCL_INCLUDE_DIRS} ${OpenCV_INCLUDE_DIRS}) +link_directories(${PCL_LIBRARY_DIRS} ${OpenCV_LIBRARY_DIRS}) +add_definitions(${PCL_DEFINITIONS}) +find_package(OpenCV 3 REQUIRED) + +message(STATUS "PCL_INCLUDE_DIRS:" ${PCL_INCLUDE_DIRS}) +message(STATUS "PCL_LIBRARY_DIRS:" ${PCL_LIBRARY_DIRS}) +message(STATUS "PCL_DEFINITIONS:" ${PCL_DEFINITIONS}) + +find_package(G2O REQUIRED) +include_directories(SYSTEM ${G2O_INCLUDE_DIR} ${G2O_INCLUDE_DIRS}) +link_directories(${G2O_LIBRARY_DIRS}) +# link_libraries(${G2O_LIBRARIES}) + +find_package(OpenMP) +if (OPENMP_FOUND) + set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") + set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +endif() + +find_library(VGICP_CUDA_FOUND NAMES fast_vgicp_cuda) +message(STATUS "VGICP_CUDA_FOUND:" ${VGICP_CUDA_FOUND}) +if(VGICP_CUDA_FOUND) + add_definitions(-DUSE_VGICP_CUDA) +endif() + +######################## +## message generation ## +######################## +add_message_files(FILES + FloorCoeffs.msg + ScanMatchingStatus.msg +) + +add_service_files(FILES + SaveMap.srv + DumpGraph.srv +) + +generate_messages(DEPENDENCIES std_msgs geometry_msgs) + +################################### +## catkin specific configuration ## +################################### +catkin_package( + INCLUDE_DIRS include + LIBRARIES hdl_graph_slam_nodelet +# CATKIN_DEPENDS pcl_ros roscpp sensor_msgs +# DEPENDS system_lib +) + +########### +## Build ## +########### +include_directories(include) +include_directories( + ${PCL_INCLUDE_DIRS} + ${catkin_INCLUDE_DIRS} + ${OpenCV_INCLUDE_DIRS} +) + +# nodelets +add_library(prefiltering_nodelet apps/prefiltering_nodelet.cpp) +target_link_libraries(prefiltering_nodelet + ${catkin_LIBRARIES} + ${PCL_LIBRARIES} + ${OpenCV_LIBRARIES} +) + + +add_library(floor_detection_nodelet apps/floor_detection_nodelet.cpp) +target_link_libraries(floor_detection_nodelet + ${catkin_LIBRARIES} + ${PCL_LIBRARIES} + ${OpenCV_LIBRARIES} +) +add_dependencies(floor_detection_nodelet ${PROJECT_NAME}_gencpp) + + +add_library(scan_matching_odometry_nodelet + apps/scan_matching_odometry_nodelet.cpp + src/hdl_graph_slam/registrations.cpp + src/hdl_graph_slam/imageProjection.cpp + src/hdl_graph_slam/featureAssociation.cpp +) +target_link_libraries(scan_matching_odometry_nodelet + ${catkin_LIBRARIES} + ${PCL_LIBRARIES} + ${OpenCV_LIBRARIES} +) +add_dependencies(scan_matching_odometry_nodelet ${PROJECT_NAME}_gencpp) + + +add_library(hdl_graph_slam_nodelet + apps/hdl_graph_slam_nodelet.cpp + src/hdl_graph_slam/graph_slam.cpp + src/hdl_graph_slam/keyframe.cpp + src/hdl_graph_slam/map_cloud_generator.cpp + src/hdl_graph_slam/registrations.cpp + src/hdl_graph_slam/information_matrix_calculator.cpp + src/g2o/robust_kernel_io.cpp +) +target_link_libraries(hdl_graph_slam_nodelet + ${catkin_LIBRARIES} + ${PCL_LIBRARIES} + ${G2O_TYPES_DATA} + ${G2O_CORE_LIBRARY} + ${G2O_STUFF_LIBRARY} + ${G2O_SOLVER_PCG} + ${G2O_SOLVER_CSPARSE} # be aware of that CSPARSE is released under LGPL + ${G2O_SOLVER_CHOLMOD} # be aware of that cholmod is released under GPL + ${G2O_TYPES_SLAM3D} + ${G2O_TYPES_SLAM3D_ADDONS} +) +add_dependencies(hdl_graph_slam_nodelet ${PROJECT_NAME}_gencpp) + +catkin_install_python( + PROGRAMS + src/${PROJECT_NAME}/bag_player.py + src/${PROJECT_NAME}/ford2bag.py + src/${PROJECT_NAME}/map2odom_publisher.py + DESTINATION ${CATKIN_PACKAGE_BIN_DESTINATION} +) + +install(FILES nodelet_plugins.xml + DESTINATION ${CATKIN_PACKAGE_SHARE_DESTINATION} +) + +install(TARGETS + prefiltering_nodelet + floor_detection_nodelet + scan_matching_odometry_nodelet + hdl_graph_slam_nodelet + LIBRARY DESTINATION ${CATKIN_PACKAGE_LIB_DESTINATION}) + +install(DIRECTORY include/ + DESTINATION ${CATKIN_GLOBAL_INCLUDE_DESTINATION} +) diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..612dcee --- /dev/null +++ b/LICENSE @@ -0,0 +1,25 @@ +BSD 2-Clause License + +Copyright (c) 2019, k.koide +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..386d3f8 --- /dev/null +++ b/README.md @@ -0,0 +1,53 @@ +#FD-SLAM +This is an open source ROS package for real-time 6DOF SLAM using a 3D LIDAR. + +It is based on hdl_graph_slam and the steps to run our system are same with hdl-graph-slam. + +##The names and functions do not change compared with hdl-graph-slam. + +We use a novel feature-based Lidar odometry for fast scan-matching, and use a modified GICP for keyframe matching. The backend in hdl-graph-slam is reused. We have tested this package with Velodyne (HDL32e, HDL64,VLP16) and Ouster64 sensors in indoor and outdoor environments. The corresponding configure launch files are provided. + +##key modifications +src/hdl_graph_slam/imageProjection.cpp and src/hdl_graph_slam/featureAssociation.cpp are added for fast scan matching + +As for the fast-gicp, key modifications are taken place in fast_gicp_impl.cpp----calculate_covariances() + + +## Nodelets +***FD-slam*** consists of four nodelets. + +- *prefiltering_nodelet* +- *scan_matching_odometry_nodelet* +- *floor_detection_nodelet* +- *hdl_graph_slam_nodelet* + +The input point cloud is first downsampled by *prefiltering_nodelet*, and then passed to the next nodelets. While *scan_matching_odometry_nodelet* estimates the sensor pose by iteratively applying a scan matching between consecutive frames (i.e., odometry estimation), *floor_detection_nodelet* detects floor planes by RANSAC. The estimated odometry and the detected floor planes are sent to *hdl_graph_slam*. To compensate the accumulated error of the scan matching, it performs loop detection and optimizes a pose graph which takes various constraints into account. + + + +## Parameters +All the configurable parameters are listed in *launch/****.launch* as ros params. + +## Services +- */hdl_graph_slam/dump* (hdl_graph_slam/DumpGraph) + - save all the internal data (point clouds, floor coeffs, odoms, and pose graph) to a directory. +- */hdl_graph_slam/save_map* (hdl_graph_slam/SaveMap) + - save the generated map as a PCD file. + +## Requirements +***hdl_graph_slam*** requires the following libraries: + +- OpenMP +- PCL +- g2o +- suitesparse + +The following ROS packages are required: + +- geodesy +- nmea_msgs +- pcl_ros +- [ndt_omp](https://github.com/koide3/ndt_omp) +- [fast_gicp](https://github.com/SMRT-AIST/fast_gicp) + + diff --git a/README_used in hdl_graph_slam.md b/README_used in hdl_graph_slam.md new file mode 100644 index 0000000..0c33b95 --- /dev/null +++ b/README_used in hdl_graph_slam.md @@ -0,0 +1,225 @@ +# hdl_graph_slam +***hdl_graph_slam*** is an open source ROS package for real-time 6DOF SLAM using a 3D LIDAR. It is based on 3D Graph SLAM with NDT scan matching-based odometry estimation and loop detection. It also supports several graph constraints, such as GPS, IMU acceleration (gravity vector), IMU orientation (magnetic sensor), and floor plane (detected in a point cloud). We have tested this package with Velodyne (HDL32e, VLP16) and RoboSense (16 channels) sensors in indoor and outdoor environments. + + + +[video](https://drive.google.com/open?id=0B9f5zFkpn4soSG96Tkt4SFFTbms) + +[![Build](https://github.com/koide3/hdl_graph_slam/actions/workflows/build.yml/badge.svg)](https://github.com/koide3/hdl_graph_slam/actions/workflows/build.yml) on melodic & noetic + +## Nodelets +***hdl_graph_slam*** consists of four nodelets. + +- *prefiltering_nodelet* +- *scan_matching_odometry_nodelet* +- *floor_detection_nodelet* +- *hdl_graph_slam_nodelet* + +The input point cloud is first downsampled by *prefiltering_nodelet*, and then passed to the next nodelets. While *scan_matching_odometry_nodelet* estimates the sensor pose by iteratively applying a scan matching between consecutive frames (i.e., odometry estimation), *floor_detection_nodelet* detects floor planes by RANSAC. The estimated odometry and the detected floor planes are sent to *hdl_graph_slam*. To compensate the accumulated error of the scan matching, it performs loop detection and optimizes a pose graph which takes various constraints into account. + + + +## Constraints (Edges) + +You can enable/disable each constraint by changing params in the launch file, and you can also change the weight (\*_stddev) and the robust kernel (\*_robust_kernel) of each constraint. + +- ***Odometry*** + +- ***Loop closure*** + +- ***GPS*** + - */gps/geopoint* (geographic_msgs/GeoPointStamped) + - */gps/navsat* (sensor_msgs/NavSatFix) + - */gpsimu_driver/nmea_sentence* (nmea_msgs/Sentence) + +hdl_graph_slam supports several GPS message types. All the supported types contain (latitude, longitude, and altitude). hdl_graph_slam converts them into [the UTM coordinate](http://wiki.ros.org/geodesy), and adds them into the graph as 3D position constraints. If altitude is set to NaN, the GPS data is treated as a 2D constrait. GeoPoint is the most basic one, which consists of only (lat, lon, alt). Although NavSatFix provides many information, we use only (lat, lon, alt) and ignore all other data. If you're using HDL32e, you can directly connect *hdl_graph_slam* with *velodyne_driver* via */gpsimu_driver/nmea_sentence*. + +- ***IMU acceleration (gravity vector)*** + - */gpsimu_driver/imu_data* (sensor_msgs/Imu) + +This constraint rotates each pose node so that the acceleration vector associated with the node becomes vertical (as the gravity vector). This is useful to compensate for accumulated tilt rotation errors of the scan matching. Since we ignore acceleration by sensor motion, you should not give a big weight for this constraint. + +- ***IMU orientation (magnetic sensor)*** + - */gpsimu_driver/imu_data* (sensor_msgs/Imu) + + If your IMU has a reliable magnetic orientation sensor, you can add orientation data to the graph as 3D rotation constraints. Note that, magnetic orientation sensors can be affected by external magnetic disturbances. In such cases, this constraint should be disabled. + +- ***Floor plane*** + - */floor_detection/floor_coeffs* (hdl_graph_slam/FloorCoeffs) + +This constraint optimizes the graph so that the floor planes (detected by RANSAC) of the pose nodes becomes the same. This is designed to compensate the accumulated rotation error of the scan matching in large flat indoor environments. + + +## Parameters +All the configurable parameters are listed in *launch/hdl_graph_slam.launch* as ros params. + +## Services +- */hdl_graph_slam/dump* (hdl_graph_slam/DumpGraph) + - save all the internal data (point clouds, floor coeffs, odoms, and pose graph) to a directory. +- */hdl_graph_slam/save_map* (hdl_graph_slam/SaveMap) + - save the generated map as a PCD file. + +## Requirements +***hdl_graph_slam*** requires the following libraries: + +- OpenMP +- PCL +- g2o +- suitesparse + +The following ROS packages are required: + +- geodesy +- nmea_msgs +- pcl_ros +- [ndt_omp](https://github.com/koide3/ndt_omp) +- [fast_gicp](https://github.com/SMRT-AIST/fast_gicp) + +```bash +# for melodic +sudo apt-get install ros-melodic-geodesy ros-melodic-pcl-ros ros-melodic-nmea-msgs ros-melodic-libg2o +cd catkin_ws/src +git clone https://github.com/koide3/ndt_omp.git -b melodic +git clone https://github.com/SMRT-AIST/fast_gicp.git --recursive +git clone https://github.com/koide3/hdl_graph_slam + +cd .. && catkin_make -DCMAKE_BUILD_TYPE=Release + +# for noetic +sudo apt-get install ros-noetic-geodesy ros-noetic-pcl-ros ros-noetic-nmea-msgs ros-noetic-libg2o + +cd catkin_ws/src +git clone https://github.com/koide3/ndt_omp.git +git clone https://github.com/SMRT-AIST/fast_gicp.git --recursive +git clone https://github.com/koide3/hdl_graph_slam + +cd .. && catkin_make -DCMAKE_BUILD_TYPE=Release +``` + +**[optional]** *bag_player.py* script requires ProgressBar2. +```bash +sudo pip install ProgressBar2 +``` + +## Example1 (Indoor) + +Bag file (recorded in a small room): + +- [hdl_501.bag.tar.gz](http://www.aisl.cs.tut.ac.jp/databases/hdl_graph_slam/hdl_501.bag.tar.gz) (raw data, 344MB) +- [hdl_501_filtered.bag.tar.gz](http://www.aisl.cs.tut.ac.jp/databases/hdl_graph_slam/hdl_501_filtered.bag.tar.gz) (downsampled data, 57MB, **Recommended!**) + +```bash +rosparam set use_sim_time true +roslaunch hdl_graph_slam hdl_graph_slam_501.launch +``` + +```bash +roscd hdl_graph_slam/rviz +rviz -d hdl_graph_slam.rviz +``` + +```bash +rosbag play --clock hdl_501_filtered.bag +``` + +We also provide bag_player.py which automatically adjusts the playback speed and processes data as fast as possible. + +```bash +rosrun hdl_graph_slam bag_player.py hdl_501_filtered.bag +``` + +You'll see a point cloud like: + + + +You can save the generated map by: +```bash +rosservice call /hdl_graph_slam/save_map "resolution: 0.05 +destination: '/full_path_directory/map.pcd'" +``` + +## Example2 (Outdoor) + +Bag file (recorded in an outdoor environment): +- [hdl_400.bag.tar.gz](http://www.aisl.cs.tut.ac.jp/databases/hdl_graph_slam/hdl_400.bag.tar.gz) (raw data, about 900MB) + +```bash +rosparam set use_sim_time true +roslaunch hdl_graph_slam hdl_graph_slam_400.launch +``` + +```bash +roscd hdl_graph_slam/rviz +rviz -d hdl_graph_slam.rviz +``` + +```bash +rosbag play --clock hdl_400.bag +``` + + + +## Example with GPS +Ford Campus Vision and Lidar Data Set [\[URL\]](http://robots.engin.umich.edu/SoftwareData/Ford) + +The following script converts the Ford Lidar Dataset to a rosbag and plays it. In this example, ***hdl_graph_slam*** utilizes the GPS data to correct the pose graph. + +```bash +cd IJRR-Dataset-2 +rosrun hdl_graph_slam ford2bag.py dataset-2.bag +rosrun hdl_graph_slam bag_player.py dataset-2.bag +``` + + + +## Use hdl_graph_slam in your system + +1. Define the transformation between your sensors (LIDAR, IMU, GPS) and base_link of your system using static_transform_publisher (see line #11, hdl_graph_slam.launch). All the sensor data will be transformed into the common base_link frame, and then fed to the SLAM algorithm. + +2. Remap the point cloud topic of ***prefiltering_nodelet***. Like: + +```bash + + ... +``` + +## Common Problems + +### Parameter tuning guide + +The mapping quality largely depends on the parameter setting. In particular, scan matching parameters have a big impact on the result. Tune the parameters accoding to the following instructions: + +- ***registration_method*** + **[updated] In short, use FAST_GICP for most cases and FAST_VGICP or NDT_OMP if the processing speed matters** This parameter allows to change the registration method to be used for odometry estimation and loop detection. Note that GICP in PCL1.7 (ROS kinetic) or earlier has a bug in the initial guess handling. **If you are on ROS kinectic or earlier, do not use GICP**. + +- ***ndt_resolution*** + This parameter decides the voxel size of NDT. Typically larger values are good for outdoor environements (0.5 - 2.0 [m] for indoor, 2.0 - 10.0 [m] for outdoor). If you chose NDT or NDT_OMP, tweak this parameter so you can obtain a good odometry estimation result. + +- ***other parameters*** + All the configurable parameters are available in the launch file. Copy a template launch file (hdl_graph_slam_501.launch for indoor, hdl_graph_slam_400.launch for outdoor) and tweak parameters in the launch file to adapt it to your application. + +## License + +This package is released under the BSD-2-Clause License. + + +Note that the cholmod solver in g2o is licensed under GPL. You may need to build g2o without cholmod dependency to avoid the GPL. + +## Related packages + +- [interactive_slam](https://github.com/koide3/interactive_slam) +- [hdl_graph_slam](https://github.com/koide3/hdl_graph_slam) +- [hdl_localization](https://github.com/koide3/hdl_localization) +- [hdl_people_tracking](https://github.com/koide3/hdl_people_tracking) + + + +## Papers +Kenji Koide, Jun Miura, and Emanuele Menegatti, A Portable 3D LIDAR-based System for Long-term and Wide-area People Behavior Measurement, Advanced Robotic Systems, 2019 [[link]](https://www.researchgate.net/publication/331283709_A_Portable_3D_LIDAR-based_System_for_Long-term_and_Wide-area_People_Behavior_Measurement). + +## Contact +Kenji Koide, k.koide@aist.go.jp, https://staff.aist.go.jp/k.koide + +Active Intelligent Systems Laboratory, Toyohashi University of Technology, Japan [\[URL\]](http://www.aisl.cs.tut.ac.jp) +Mobile Robotics Research Team, National Institute of Advanced Industrial Science and Technology (AIST), Japan [\[URL\]](https://unit.aist.go.jp/rirc/en/team/smart_mobility.html) diff --git a/apps/floor_detection_nodelet.cpp b/apps/floor_detection_nodelet.cpp new file mode 100644 index 0000000..ad714d9 --- /dev/null +++ b/apps/floor_detection_nodelet.cpp @@ -0,0 +1,269 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace hdl_graph_slam { + +class FloorDetectionNodelet : public nodelet::Nodelet { +public: + typedef pcl::PointXYZI PointT; + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + + FloorDetectionNodelet() {} + virtual ~FloorDetectionNodelet() {} + + virtual void onInit() { + NODELET_DEBUG("initializing floor_detection_nodelet..."); + nh = getNodeHandle(); + private_nh = getPrivateNodeHandle(); + + initialize_params(); + + points_sub = nh.subscribe("/filtered_points", 256, &FloorDetectionNodelet::cloud_callback, this); + floor_pub = nh.advertise("/floor_detection/floor_coeffs", 32); + + read_until_pub = nh.advertise("/floor_detection/read_until", 32); + floor_filtered_pub = nh.advertise("/floor_detection/floor_filtered_points", 32); + floor_points_pub = nh.advertise("/floor_detection/floor_points", 32); + } + +private: + /** + * @brief initialize parameters + */ + void initialize_params() { + tilt_deg = private_nh.param("tilt_deg", 0.0); // approximate sensor tilt angle [deg] + sensor_height = private_nh.param("sensor_height", 2.0); // approximate sensor height [m] + height_clip_range = private_nh.param("height_clip_range", 1.0); // points with heights in [sensor_height - height_clip_range, sensor_height + height_clip_range] will be used for floor detection + floor_pts_thresh = private_nh.param("floor_pts_thresh", 512); // minimum number of support points of RANSAC to accept a detected floor plane + floor_normal_thresh = private_nh.param("floor_normal_thresh", 10.0); // verticality check thresold for the detected floor plane [deg] + use_normal_filtering = private_nh.param("use_normal_filtering", true); // if true, points with "non-"vertical normals will be filtered before RANSAC + normal_filter_thresh = private_nh.param("normal_filter_thresh", 20.0); // "non-"verticality check threshold [deg] + + points_topic = private_nh.param("points_topic", "/velodyne_points"); + } + + /** + * @brief callback for point clouds + * @param cloud_msg point cloud msg + */ + void cloud_callback(const sensor_msgs::PointCloud2ConstPtr& cloud_msg) { + pcl::PointCloud::Ptr cloud(new pcl::PointCloud()); + pcl::fromROSMsg(*cloud_msg, *cloud); + + if(cloud->empty()) { + return; + } + + // floor detection + boost::optional floor = detect(cloud); + + // publish the detected floor coefficients + hdl_graph_slam::FloorCoeffs coeffs; + coeffs.header = cloud_msg->header; + if(floor) { + coeffs.coeffs.resize(4); + for(int i = 0; i < 4; i++) { + coeffs.coeffs[i] = (*floor)[i]; + } + } + + floor_pub.publish(coeffs); + + // for offline estimation + std_msgs::HeaderPtr read_until(new std_msgs::Header()); + read_until->frame_id = points_topic; + read_until->stamp = cloud_msg->header.stamp + ros::Duration(1, 0); + read_until_pub.publish(read_until); + + read_until->frame_id = "/filtered_points"; + read_until_pub.publish(read_until); + } + + /** + * @brief detect the floor plane from a point cloud + * @param cloud input cloud + * @return detected floor plane coefficients + */ + boost::optional detect(const pcl::PointCloud::Ptr& cloud) const { + // compensate the tilt rotation + Eigen::Matrix4f tilt_matrix = Eigen::Matrix4f::Identity(); + tilt_matrix.topLeftCorner(3, 3) = Eigen::AngleAxisf(tilt_deg * M_PI / 180.0f, Eigen::Vector3f::UnitY()).toRotationMatrix(); + + // filtering before RANSAC (height and normal filtering) + pcl::PointCloud::Ptr filtered(new pcl::PointCloud); + pcl::transformPointCloud(*cloud, *filtered, tilt_matrix); + filtered = plane_clip(filtered, Eigen::Vector4f(0.0f, 0.0f, 1.0f, sensor_height + height_clip_range), false); + filtered = plane_clip(filtered, Eigen::Vector4f(0.0f, 0.0f, 1.0f, sensor_height - height_clip_range), true); + + if(use_normal_filtering) { + filtered = normal_filtering(filtered); + } + + pcl::transformPointCloud(*filtered, *filtered, static_cast(tilt_matrix.inverse())); + + if(floor_filtered_pub.getNumSubscribers()) { + filtered->header = cloud->header; + floor_filtered_pub.publish(*filtered); + } + + // too few points for RANSAC + if(filtered->size() < floor_pts_thresh) { + return boost::none; + } + + // RANSAC + pcl::SampleConsensusModelPlane::Ptr model_p(new pcl::SampleConsensusModelPlane(filtered)); + pcl::RandomSampleConsensus ransac(model_p); + ransac.setDistanceThreshold(0.1); + ransac.computeModel(); + + pcl::PointIndices::Ptr inliers(new pcl::PointIndices); + ransac.getInliers(inliers->indices); + + // too few inliers + if(inliers->indices.size() < floor_pts_thresh) { + return boost::none; + } + + // verticality check of the detected floor's normal + Eigen::Vector4f reference = tilt_matrix.inverse() * Eigen::Vector4f::UnitZ(); + + Eigen::VectorXf coeffs; + ransac.getModelCoefficients(coeffs); + + double dot = coeffs.head<3>().dot(reference.head<3>()); + if(std::abs(dot) < std::cos(floor_normal_thresh * M_PI / 180.0)) { + // the normal is not vertical + return boost::none; + } + + // make the normal upward + if(coeffs.head<3>().dot(Eigen::Vector3f::UnitZ()) < 0.0f) { + coeffs *= -1.0f; + } + + if(floor_points_pub.getNumSubscribers()) { + pcl::PointCloud::Ptr inlier_cloud(new pcl::PointCloud); + pcl::ExtractIndices extract; + extract.setInputCloud(filtered); + extract.setIndices(inliers); + extract.filter(*inlier_cloud); + inlier_cloud->header = cloud->header; + + floor_points_pub.publish(*inlier_cloud); + } + + return Eigen::Vector4f(coeffs); + } + + /** + * @brief plane_clip + * @param src_cloud + * @param plane + * @param negative + * @return + */ + pcl::PointCloud::Ptr plane_clip(const pcl::PointCloud::Ptr& src_cloud, const Eigen::Vector4f& plane, bool negative) const { + pcl::PlaneClipper3D clipper(plane); + pcl::PointIndices::Ptr indices(new pcl::PointIndices); + + clipper.clipPointCloud3D(*src_cloud, indices->indices); + + pcl::PointCloud::Ptr dst_cloud(new pcl::PointCloud); + + pcl::ExtractIndices extract; + extract.setInputCloud(src_cloud); + extract.setIndices(indices); + extract.setNegative(negative); + extract.filter(*dst_cloud); + + return dst_cloud; + } + + /** + * @brief filter points with non-vertical normals + * @param cloud input cloud + * @return filtered cloud + */ + pcl::PointCloud::Ptr normal_filtering(const pcl::PointCloud::Ptr& cloud) const { + pcl::NormalEstimation ne; + ne.setInputCloud(cloud); + + pcl::search::KdTree::Ptr tree(new pcl::search::KdTree); + ne.setSearchMethod(tree); + + pcl::PointCloud::Ptr normals(new pcl::PointCloud); + ne.setKSearch(10); + ne.setViewPoint(0.0f, 0.0f, sensor_height); + ne.compute(*normals); + + pcl::PointCloud::Ptr filtered(new pcl::PointCloud); + filtered->reserve(cloud->size()); + + for(int i = 0; i < cloud->size(); i++) { + float dot = normals->at(i).getNormalVector3fMap().normalized().dot(Eigen::Vector3f::UnitZ()); + if(std::abs(dot) > std::cos(normal_filter_thresh * M_PI / 180.0)) { + filtered->push_back(cloud->at(i)); + } + } + + filtered->width = filtered->size(); + filtered->height = 1; + filtered->is_dense = false; + + return filtered; + } + +private: + ros::NodeHandle nh; + ros::NodeHandle private_nh; + + // ROS topics + ros::Subscriber points_sub; + + ros::Publisher floor_pub; + ros::Publisher floor_points_pub; + ros::Publisher floor_filtered_pub; + + std::string points_topic; + ros::Publisher read_until_pub; + + // floor detection parameters + // see initialize_params() for the details + double tilt_deg; + double sensor_height; + double height_clip_range; + + int floor_pts_thresh; + double floor_normal_thresh; + + bool use_normal_filtering; + double normal_filter_thresh; +}; + +} // namespace hdl_graph_slam + +PLUGINLIB_EXPORT_CLASS(hdl_graph_slam::FloorDetectionNodelet, nodelet::Nodelet) diff --git a/apps/hdl_graph_slam_nodelet.cpp b/apps/hdl_graph_slam_nodelet.cpp new file mode 100644 index 0000000..23c9af2 --- /dev/null +++ b/apps/hdl_graph_slam_nodelet.cpp @@ -0,0 +1,1056 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +using namespace std; +namespace hdl_graph_slam { + +class HdlGraphSlamNodelet : public nodelet::Nodelet { +public: + typedef pcl::PointXYZI PointT; + typedef message_filters::sync_policies::ApproximateTime ApproxSyncPolicy; + + HdlGraphSlamNodelet() {} + virtual ~HdlGraphSlamNodelet() {} + + virtual void onInit() { + nh = getNodeHandle(); + mt_nh = getMTNodeHandle(); + private_nh = getPrivateNodeHandle(); + + // init parameters + map_frame_id = private_nh.param("map_frame_id", "map"); + odom_frame_id = private_nh.param("odom_frame_id", "odom"); + map_cloud_resolution = private_nh.param("map_cloud_resolution", 0.05); + trans_odom2map.setIdentity(); + + max_keyframes_per_update = private_nh.param("max_keyframes_per_update", 10); + + // + anchor_node = nullptr; + anchor_edge = nullptr; + floor_plane_node = nullptr; + graph_slam.reset(new GraphSLAM(private_nh.param("g2o_solver_type", "lm_var"))); + keyframe_updater.reset(new KeyframeUpdater(private_nh)); + loop_detector.reset(new LoopDetector(private_nh)); + map_cloud_generator.reset(new MapCloudGenerator()); + inf_calclator.reset(new InformationMatrixCalculator(private_nh)); + nmea_parser.reset(new NmeaSentenceParser()); + + gps_time_offset = private_nh.param("gps_time_offset", 0.0); + gps_edge_stddev_xy = private_nh.param("gps_edge_stddev_xy", 10000.0); + gps_edge_stddev_z = private_nh.param("gps_edge_stddev_z", 10.0); + floor_edge_stddev = private_nh.param("floor_edge_stddev", 10.0); + seq_ = private_nh.param("Sequence", "seq"); + std::cout<<"seq: "<("use_prev_test_", false); + + imu_time_offset = private_nh.param("imu_time_offset", 0.0); + enable_imu_orientation = private_nh.param("enable_imu_orientation", false); + enable_imu_acceleration = private_nh.param("enable_imu_acceleration", false); + imu_orientation_edge_stddev = private_nh.param("imu_orientation_edge_stddev", 0.1); + imu_acceleration_edge_stddev = private_nh.param("imu_acceleration_edge_stddev", 3.0); + + points_topic = private_nh.param("points_topic", "/velodyne_points"); + + // subscribers + odom_sub.reset(new message_filters::Subscriber(mt_nh, "/odom", 256)); + cloud_sub.reset(new message_filters::Subscriber(mt_nh, "/filtered_points", 32)); + sync.reset(new message_filters::Synchronizer(ApproxSyncPolicy(32), *odom_sub, *cloud_sub)); + sync->registerCallback(boost::bind(&HdlGraphSlamNodelet::cloud_callback, this, _1, _2)); + imu_sub = nh.subscribe("/gpsimu_driver/imu_data", 1024, &HdlGraphSlamNodelet::imu_callback, this); + floor_sub = nh.subscribe("/floor_detection/floor_coeffs", 1024, &HdlGraphSlamNodelet::floor_coeffs_callback, this); + + if(private_nh.param("enable_gps", true)) { + gps_sub = mt_nh.subscribe("/gps/geopoint", 1024, &HdlGraphSlamNodelet::gps_callback, this); + nmea_sub = mt_nh.subscribe("/gpsimu_driver/nmea_sentence", 1024, &HdlGraphSlamNodelet::nmea_callback, this); + navsat_sub = mt_nh.subscribe("/gps/navsat", 1024, &HdlGraphSlamNodelet::navsat_callback, this); + } + + // publishers + markers_pub = mt_nh.advertise("/hdl_graph_slam/markers", 16); + odom2map_pub = mt_nh.advertise("/hdl_graph_slam/odom2pub", 16); + map_points_pub = mt_nh.advertise("/hdl_graph_slam/map_points", 1, true); + read_until_pub = mt_nh.advertise("/hdl_graph_slam/read_until", 32); + + dump_service_server = mt_nh.advertiseService("/hdl_graph_slam/dump", &HdlGraphSlamNodelet::dump_service, this); + save_map_service_server = mt_nh.advertiseService("/hdl_graph_slam/save_map", &HdlGraphSlamNodelet::save_map_service, this); + + graph_updated = false; + double graph_update_interval = private_nh.param("graph_update_interval", 3.0); + double map_cloud_update_interval = private_nh.param("map_cloud_update_interval", 10.0); + optimization_timer = mt_nh.createWallTimer(ros::WallDuration(graph_update_interval), &HdlGraphSlamNodelet::optimization_timer_callback, this); + map_publish_timer = mt_nh.createWallTimer(ros::WallDuration(map_cloud_update_interval), &HdlGraphSlamNodelet::map_points_publish_timer_callback, this); + } + +private: + /** + * @brief received point clouds are pushed to #keyframe_queue + * @param odom_msg + * @param cloud_msg + */ + void cloud_callback(const nav_msgs::OdometryConstPtr& odom_msg, const sensor_msgs::PointCloud2::ConstPtr& cloud_msg) { + const ros::Time& stamp = cloud_msg->header.stamp; + Eigen::Isometry3d odom = odom2isometry(odom_msg); + + pcl::PointCloud::Ptr cloud(new pcl::PointCloud()); + pcl::fromROSMsg(*cloud_msg, *cloud); + if(base_frame_id.empty()) { + base_frame_id = cloud_msg->header.frame_id; + } + + if(!keyframe_updater->update(odom)) { + std::lock_guard lock(keyframe_queue_mutex); + if(keyframe_queue.empty()) { + std_msgs::Header read_until; + read_until.stamp = stamp + ros::Duration(10, 0); + read_until.frame_id = points_topic; + read_until_pub.publish(read_until); + read_until.frame_id = "/filtered_points"; + read_until_pub.publish(read_until); + } + + return; + } + + double accum_d = keyframe_updater->get_accum_distance(); + + pcl::Filter::Ptr keyframe_downsample_filter; + double res_ = 0.5; + auto keyframe_voxelgrid = new pcl::VoxelGrid(); + keyframe_voxelgrid->setLeafSize(res_, res_, res_); + keyframe_downsample_filter.reset(keyframe_voxelgrid); + + pcl::PointCloud::Ptr filtered(new pcl::PointCloud()); + keyframe_downsample_filter->setInputCloud(cloud); + keyframe_downsample_filter->filter(*filtered); + KeyFrame::Ptr keyframe(new KeyFrame(stamp, odom, accum_d, filtered)); + + std::lock_guard lock(keyframe_queue_mutex); + keyframe_queue.push_back(keyframe); + } + + /** + * @brief this method adds all the keyframes in #keyframe_queue to the pose graph (odometry edges) + * @return if true, at least one keyframe was added to the pose graph + */ + bool flush_keyframe_queue() { + std::lock_guard lock(keyframe_queue_mutex); + + if(keyframe_queue.empty()) { + return false; + } + + trans_odom2map_mutex.lock(); + Eigen::Isometry3d odom2map(trans_odom2map.cast()); + trans_odom2map_mutex.unlock(); + + int num_processed = 0; + for(int i = 0; i < std::min(keyframe_queue.size(), max_keyframes_per_update); i++) { + num_processed = i; + + const auto& keyframe = keyframe_queue[i]; + // new_keyframes will be tested later for loop closure + new_keyframes.push_back(keyframe); + + // add pose node + Eigen::Isometry3d odom = odom2map * keyframe->odom; + keyframe->node = graph_slam->add_se3_node(odom); + keyframe_hash[keyframe->stamp] = keyframe; + + // fix the first node + if(keyframes.empty() && new_keyframes.size() == 1) { + if(private_nh.param("fix_first_node", false)) { + Eigen::MatrixXd inf = Eigen::MatrixXd::Identity(6, 6); + std::stringstream sst(private_nh.param("fix_first_node_stddev", "1 1 1 1 1 1")); + for(int i = 0; i < 6; i++) { + double stddev = 1.0; + sst >> stddev; + inf(i, i) = 1.0 / stddev; + } + + anchor_node = graph_slam->add_se3_node(Eigen::Isometry3d::Identity()); + anchor_node->setFixed(true); + anchor_edge = graph_slam->add_se3_edge(anchor_node, keyframe->node, Eigen::Isometry3d::Identity(), inf); + } + } + + if(i == 0 && keyframes.empty()) { + continue; + } + + // add edge between consecutive keyframes + const auto& prev_keyframe = i == 0 ? keyframes.back() : keyframe_queue[i - 1]; + + Eigen::Isometry3d relative_pose = keyframe->odom.inverse() * prev_keyframe->odom; + Eigen::MatrixXd information = inf_calclator->calc_information_matrix(keyframe->cloud, prev_keyframe->cloud, relative_pose); + auto edge = graph_slam->add_se3_edge(keyframe->node, prev_keyframe->node, relative_pose, information); + graph_slam->add_robust_kernel(edge, private_nh.param("odometry_edge_robust_kernel", "NONE"), private_nh.param("odometry_edge_robust_kernel_size", 1.0)); + } + + std_msgs::Header read_until; + read_until.stamp = keyframe_queue[num_processed]->stamp + ros::Duration(10, 0); + read_until.frame_id = points_topic; + read_until_pub.publish(read_until); + read_until.frame_id = "/filtered_points"; + read_until_pub.publish(read_until); + + keyframe_queue.erase(keyframe_queue.begin(), keyframe_queue.begin() + num_processed + 1); + return true; + } + + void nmea_callback(const nmea_msgs::SentenceConstPtr& nmea_msg) { + GPRMC grmc = nmea_parser->parse(nmea_msg->sentence); + + if(grmc.status != 'A') { + return; + } + + geographic_msgs::GeoPointStampedPtr gps_msg(new geographic_msgs::GeoPointStamped()); + gps_msg->header = nmea_msg->header; + gps_msg->position.latitude = grmc.latitude; + gps_msg->position.longitude = grmc.longitude; + gps_msg->position.altitude = NAN; + + gps_callback(gps_msg); + } + + void navsat_callback(const sensor_msgs::NavSatFixConstPtr& navsat_msg) { + geographic_msgs::GeoPointStampedPtr gps_msg(new geographic_msgs::GeoPointStamped()); + gps_msg->header = navsat_msg->header; + gps_msg->position.latitude = navsat_msg->latitude; + gps_msg->position.longitude = navsat_msg->longitude; + gps_msg->position.altitude = navsat_msg->altitude; + gps_callback(gps_msg); + } + + /** + * @brief received gps data is added to #gps_queue + * @param gps_msg + */ + void gps_callback(const geographic_msgs::GeoPointStampedPtr& gps_msg) { + std::lock_guard lock(gps_queue_mutex); + gps_msg->header.stamp += ros::Duration(gps_time_offset); + gps_queue.push_back(gps_msg); + } + + /** + * @brief + * @return + */ + bool flush_gps_queue() { + std::lock_guard lock(gps_queue_mutex); + + if(keyframes.empty() || gps_queue.empty()) { + return false; + } + + bool updated = false; + auto gps_cursor = gps_queue.begin(); + + for(auto& keyframe : keyframes) { + if(keyframe->stamp > gps_queue.back()->header.stamp) { + break; + } + + if(keyframe->stamp < (*gps_cursor)->header.stamp || keyframe->utm_coord) { + continue; + } + + // find the gps data which is closest to the keyframe + auto closest_gps = gps_cursor; + for(auto gps = gps_cursor; gps != gps_queue.end(); gps++) { + auto dt = ((*closest_gps)->header.stamp - keyframe->stamp).toSec(); + auto dt2 = ((*gps)->header.stamp - keyframe->stamp).toSec(); + if(std::abs(dt) < std::abs(dt2)) { + break; + } + + closest_gps = gps; + } + + // if the time residual between the gps and keyframe is too large, skip it + gps_cursor = closest_gps; + if(0.2 < std::abs(((*closest_gps)->header.stamp - keyframe->stamp).toSec())) { + continue; + } + + // convert (latitude, longitude, altitude) -> (easting, northing, altitude) in UTM coordinate + geodesy::UTMPoint utm; + geodesy::fromMsg((*closest_gps)->position, utm); + Eigen::Vector3d xyz(utm.easting, utm.northing, utm.altitude); + + // the first gps data position will be the origin of the map + if(!zero_utm) { + zero_utm = xyz; + } + xyz -= (*zero_utm); + + keyframe->utm_coord = xyz; + + g2o::OptimizableGraph::Edge* edge; + if(std::isnan(xyz.z())) { + Eigen::Matrix2d information_matrix = Eigen::Matrix2d::Identity() / gps_edge_stddev_xy; + edge = graph_slam->add_se3_prior_xy_edge(keyframe->node, xyz.head<2>(), information_matrix); + } else { + Eigen::Matrix3d information_matrix = Eigen::Matrix3d::Identity(); + information_matrix.block<2, 2>(0, 0) /= gps_edge_stddev_xy; + information_matrix(2, 2) /= gps_edge_stddev_z; + edge = graph_slam->add_se3_prior_xyz_edge(keyframe->node, xyz, information_matrix); + } + graph_slam->add_robust_kernel(edge, private_nh.param("gps_edge_robust_kernel", "NONE"), private_nh.param("gps_edge_robust_kernel_size", 1.0)); + + updated = true; + } + + auto remove_loc = std::upper_bound(gps_queue.begin(), gps_queue.end(), keyframes.back()->stamp, [=](const ros::Time& stamp, const geographic_msgs::GeoPointStampedConstPtr& geopoint) { return stamp < geopoint->header.stamp; }); + gps_queue.erase(gps_queue.begin(), remove_loc); + return updated; + } + + void imu_callback(const sensor_msgs::ImuPtr& imu_msg) { + if(!enable_imu_orientation && !enable_imu_acceleration) { + return; + } + + std::lock_guard lock(imu_queue_mutex); + imu_msg->header.stamp += ros::Duration(imu_time_offset); + imu_queue.push_back(imu_msg); + } + + bool flush_imu_queue() { + std::lock_guard lock(imu_queue_mutex); + if(keyframes.empty() || imu_queue.empty() || base_frame_id.empty()) { + return false; + } + + bool updated = false; + auto imu_cursor = imu_queue.begin(); + + for(auto& keyframe : keyframes) { + if(keyframe->stamp > imu_queue.back()->header.stamp) { + break; + } + + if(keyframe->stamp < (*imu_cursor)->header.stamp || keyframe->acceleration) { + continue; + } + + // find imu data which is closest to the keyframe + auto closest_imu = imu_cursor; + for(auto imu = imu_cursor; imu != imu_queue.end(); imu++) { + auto dt = ((*closest_imu)->header.stamp - keyframe->stamp).toSec(); + auto dt2 = ((*imu)->header.stamp - keyframe->stamp).toSec(); + if(std::abs(dt) < std::abs(dt2)) { + break; + } + + closest_imu = imu; + } + + imu_cursor = closest_imu; + if(0.2 < std::abs(((*closest_imu)->header.stamp - keyframe->stamp).toSec())) { + continue; + } + + const auto& imu_ori = (*closest_imu)->orientation; + const auto& imu_acc = (*closest_imu)->linear_acceleration; + + geometry_msgs::Vector3Stamped acc_imu; + geometry_msgs::Vector3Stamped acc_base; + geometry_msgs::QuaternionStamped quat_imu; + geometry_msgs::QuaternionStamped quat_base; + + quat_imu.header.frame_id = acc_imu.header.frame_id = (*closest_imu)->header.frame_id; + quat_imu.header.stamp = acc_imu.header.stamp = ros::Time(0); + acc_imu.vector = (*closest_imu)->linear_acceleration; + quat_imu.quaternion = (*closest_imu)->orientation; + + try { + tf_listener.transformVector(base_frame_id, acc_imu, acc_base); + tf_listener.transformQuaternion(base_frame_id, quat_imu, quat_base); + } catch(std::exception& e) { + std::cerr << "failed to find transform!!" << std::endl; + return false; + } + + keyframe->acceleration = Eigen::Vector3d(acc_base.vector.x, acc_base.vector.y, acc_base.vector.z); + keyframe->orientation = Eigen::Quaterniond(quat_base.quaternion.w, quat_base.quaternion.x, quat_base.quaternion.y, quat_base.quaternion.z); + keyframe->orientation = keyframe->orientation; + if(keyframe->orientation->w() < 0.0) { + keyframe->orientation->coeffs() = -keyframe->orientation->coeffs(); + } + + if(enable_imu_orientation) { + Eigen::MatrixXd info = Eigen::MatrixXd::Identity(3, 3) / imu_orientation_edge_stddev; + auto edge = graph_slam->add_se3_prior_quat_edge(keyframe->node, *keyframe->orientation, info); + graph_slam->add_robust_kernel(edge, private_nh.param("imu_orientation_edge_robust_kernel", "NONE"), private_nh.param("imu_orientation_edge_robust_kernel_size", 1.0)); + } + + if(enable_imu_acceleration) { + Eigen::MatrixXd info = Eigen::MatrixXd::Identity(3, 3) / imu_acceleration_edge_stddev; + g2o::OptimizableGraph::Edge* edge = graph_slam->add_se3_prior_vec_edge(keyframe->node, -Eigen::Vector3d::UnitZ(), *keyframe->acceleration, info); + graph_slam->add_robust_kernel(edge, private_nh.param("imu_acceleration_edge_robust_kernel", "NONE"), private_nh.param("imu_acceleration_edge_robust_kernel_size", 1.0)); + } + updated = true; + } + + auto remove_loc = std::upper_bound(imu_queue.begin(), imu_queue.end(), keyframes.back()->stamp, [=](const ros::Time& stamp, const sensor_msgs::ImuConstPtr& imu) { return stamp < imu->header.stamp; }); + imu_queue.erase(imu_queue.begin(), remove_loc); + + return updated; + } + + /** + * @brief received floor coefficients are added to #floor_coeffs_queue + * @param floor_coeffs_msg + */ + void floor_coeffs_callback(const hdl_graph_slam::FloorCoeffsConstPtr& floor_coeffs_msg) { + if(floor_coeffs_msg->coeffs.empty()) { + return; + } + + std::lock_guard lock(floor_coeffs_queue_mutex); + floor_coeffs_queue.push_back(floor_coeffs_msg); + } + + /** + * @brief this methods associates floor coefficients messages with registered keyframes, and then adds the associated coeffs to the pose graph + * @return if true, at least one floor plane edge is added to the pose graph + */ + bool flush_floor_queue() { + std::lock_guard lock(floor_coeffs_queue_mutex); + + if(keyframes.empty()) { + return false; + } + + const auto& latest_keyframe_stamp = keyframes.back()->stamp; + + bool updated = false; + for(const auto& floor_coeffs : floor_coeffs_queue) { + if(floor_coeffs->header.stamp > latest_keyframe_stamp) { + break; + } + + auto found = keyframe_hash.find(floor_coeffs->header.stamp); + if(found == keyframe_hash.end()) { + continue; + } + + if(!floor_plane_node) { + floor_plane_node = graph_slam->add_plane_node(Eigen::Vector4d(0.0, 0.0, 1.0, 0.0)); + floor_plane_node->setFixed(true); + } + + const auto& keyframe = found->second; + + Eigen::Vector4d coeffs(floor_coeffs->coeffs[0], floor_coeffs->coeffs[1], floor_coeffs->coeffs[2], floor_coeffs->coeffs[3]); + Eigen::Matrix3d information = Eigen::Matrix3d::Identity() * (1.0 / floor_edge_stddev); + auto edge = graph_slam->add_se3_plane_edge(keyframe->node, floor_plane_node, coeffs, information); + graph_slam->add_robust_kernel(edge, private_nh.param("floor_edge_robust_kernel", "NONE"), private_nh.param("floor_edge_robust_kernel_size", 1.0)); + + keyframe->floor_coeffs = coeffs; + + updated = true; + } + + auto remove_loc = std::upper_bound(floor_coeffs_queue.begin(), floor_coeffs_queue.end(), latest_keyframe_stamp, [=](const ros::Time& stamp, const hdl_graph_slam::FloorCoeffsConstPtr& coeffs) { return stamp < coeffs->header.stamp; }); + floor_coeffs_queue.erase(floor_coeffs_queue.begin(), remove_loc); + return updated; + } + + /** + * @brief generate map point cloud and publish it + * @param event + */ + void map_points_publish_timer_callback(const ros::WallTimerEvent& event) { + if(!map_points_pub.getNumSubscribers() || !graph_updated) { + return; + } + + std::vector snapshot; + + keyframes_snapshot_mutex.lock(); + snapshot = keyframes_snapshot; + keyframes_snapshot_mutex.unlock(); + + auto cloud = map_cloud_generator->generate(snapshot, map_cloud_resolution); + if(!cloud) { + return; + } + + cloud->header.frame_id = map_frame_id; + cloud->header.stamp = snapshot.back()->cloud->header.stamp; + + sensor_msgs::PointCloud2Ptr cloud_msg(new sensor_msgs::PointCloud2()); + pcl::toROSMsg(*cloud, *cloud_msg); + + map_points_pub.publish(cloud_msg); + } + + /** + * @brief this methods adds all the data in the queues to the pose graph, and then optimizes the pose graph + * @param event + */ + void optimization_timer_callback(const ros::WallTimerEvent& event) { + std::lock_guard lock(main_thread_mutex); + + // add keyframes and floor coeffs in the queues to the pose graph + bool keyframe_updated = flush_keyframe_queue(); + + if(!keyframe_updated) { + std_msgs::Header read_until; + read_until.stamp = ros::Time::now() + ros::Duration(30, 0); + read_until.frame_id = points_topic; + read_until_pub.publish(read_until); + read_until.frame_id = "/filtered_points"; + read_until_pub.publish(read_until); + } + + + if(!keyframe_updated & !flush_floor_queue()) { + return; + } + + // loop detection + bool loop_detect_ = true; + if(loop_detect_) + { + std::vector loops = loop_detector->detect(keyframes, new_keyframes, *graph_slam); + std::cout<<"loops size: "<relative_pose.cast()); + Eigen::MatrixXd information_matrix = inf_calclator->calc_information_matrix(loop->key1->cloud, loop->key2->cloud, relpose); + auto edge = graph_slam->add_se3_edge(loop->key1->node, loop->key2->node, relpose, information_matrix); + graph_slam->add_robust_kernel(edge, private_nh.param("loop_closure_edge_robust_kernel", "NONE"), private_nh.param("loop_closure_edge_robust_kernel_size", 1.0)); + } + + } + + std::copy(new_keyframes.begin(), new_keyframes.end(), std::back_inserter(keyframes)); + new_keyframes.clear(); + + // move the first node anchor position to the current estimate of the first node pose + // so the first node moves freely while trying to stay around the origin + if(anchor_node && private_nh.param("fix_first_node_adaptive", true)) { + Eigen::Isometry3d anchor_target = static_cast(anchor_edge->vertices()[1])->estimate(); + anchor_node->setEstimate(anchor_target); + } + + // optimize the pose graph + bool use_g2o = true; + if(use_g2o) + { + int num_iterations = private_nh.param("g2o_solver_num_iterations", 1024); + graph_slam->optimize(num_iterations); + + } + + + //TO save the trajectory. + cout<<"begin save trajectory"<node->estimate(); + keyframe_M = keyframe_T.matrix(); + R = keyframe_M.block<3,3>(0,0); + q = R; + traj<stamp<<" "<get_accum_distance()<node->estimate() * keyframe->odom.inverse(); + trans_odom2map_mutex.lock(); + trans_odom2map = trans.matrix().cast(); + trans_odom2map_mutex.unlock(); + + std::vector snapshot(keyframes.size()); + std::transform(keyframes.begin(), keyframes.end(), snapshot.begin(), [=](const KeyFrame::Ptr& k) { return std::make_shared(k); }); + + keyframes_snapshot_mutex.lock(); + keyframes_snapshot.swap(snapshot); + keyframes_snapshot_mutex.unlock(); + graph_updated = true; + + if(odom2map_pub.getNumSubscribers()) { + geometry_msgs::TransformStamped ts = matrix2transform(keyframe->stamp, trans.matrix().cast(), map_frame_id, odom_frame_id); + odom2map_pub.publish(ts); + } + + if(markers_pub.getNumSubscribers()) { + auto markers = create_marker_array(ros::Time::now()); + markers_pub.publish(markers); + } + } + + /** + * @brief create visualization marker + * @param stamp + * @return + */ + visualization_msgs::MarkerArray create_marker_array(const ros::Time& stamp) const { + visualization_msgs::MarkerArray markers; + markers.markers.resize(4); + + // node markers + visualization_msgs::Marker& traj_marker = markers.markers[0]; + traj_marker.header.frame_id = "map"; + traj_marker.header.stamp = stamp; + traj_marker.ns = "nodes"; + traj_marker.id = 0; + traj_marker.type = visualization_msgs::Marker::SPHERE_LIST; + + traj_marker.pose.orientation.w = 1.0; + traj_marker.scale.x = traj_marker.scale.y = traj_marker.scale.z = 2.0;//1.0 + + visualization_msgs::Marker& imu_marker = markers.markers[1]; + imu_marker.header = traj_marker.header; + imu_marker.ns = "imu"; + imu_marker.id = 1; + imu_marker.type = visualization_msgs::Marker::SPHERE_LIST; + + imu_marker.pose.orientation.w = 1.0; + imu_marker.scale.x = imu_marker.scale.y = imu_marker.scale.z = 0.75; + + traj_marker.points.resize(keyframes.size()); + traj_marker.colors.resize(keyframes.size()); + for(int i = 0; i < keyframes.size(); i++) { + Eigen::Vector3d pos = keyframes[i]->node->estimate().translation(); + traj_marker.points[i].x = pos.x(); + traj_marker.points[i].y = pos.y(); + traj_marker.points[i].z = pos.z(); + + double p = static_cast(i) / keyframes.size(); + traj_marker.colors[i].r = 1.0 - p; + traj_marker.colors[i].g = p; + traj_marker.colors[i].b = 0.0; + traj_marker.colors[i].a = 1.0; + + if(keyframes[i]->acceleration) { + Eigen::Vector3d pos = keyframes[i]->node->estimate().translation(); + geometry_msgs::Point point; + point.x = pos.x(); + point.y = pos.y(); + point.z = pos.z(); + + std_msgs::ColorRGBA color; + color.r = 0.0; + color.g = 0.0; + color.b = 1.0; + color.a = 0.1; + + imu_marker.points.push_back(point); + imu_marker.colors.push_back(color); + } + } + + // edge markers + visualization_msgs::Marker& edge_marker = markers.markers[2]; + edge_marker.header.frame_id = "map"; + edge_marker.header.stamp = stamp; + edge_marker.ns = "edges"; + edge_marker.id = 2; + edge_marker.type = visualization_msgs::Marker::LINE_LIST; + + edge_marker.pose.orientation.w = 1.0; + edge_marker.scale.x = 2;//0.05 + + edge_marker.points.resize(graph_slam->graph->edges().size() * 2); + edge_marker.colors.resize(graph_slam->graph->edges().size() * 2); + + auto edge_itr = graph_slam->graph->edges().begin(); + for(int i = 0; edge_itr != graph_slam->graph->edges().end(); edge_itr++, i++) { + g2o::HyperGraph::Edge* edge = *edge_itr; + g2o::EdgeSE3* edge_se3 = dynamic_cast(edge); + if(edge_se3) { + g2o::VertexSE3* v1 = dynamic_cast(edge_se3->vertices()[0]); + g2o::VertexSE3* v2 = dynamic_cast(edge_se3->vertices()[1]); + Eigen::Vector3d pt1 = v1->estimate().translation(); + Eigen::Vector3d pt2 = v2->estimate().translation(); + + edge_marker.points[i * 2].x = pt1.x(); + edge_marker.points[i * 2].y = pt1.y(); + edge_marker.points[i * 2].z = pt1.z(); + edge_marker.points[i * 2 + 1].x = pt2.x(); + edge_marker.points[i * 2 + 1].y = pt2.y(); + edge_marker.points[i * 2 + 1].z = pt2.z(); + + double p1 = static_cast(v1->id()) / graph_slam->graph->vertices().size(); + double p2 = static_cast(v2->id()) / graph_slam->graph->vertices().size(); + + edge_marker.colors[i * 2].r = 1.0 - p1; + edge_marker.colors[i * 2].g = p1; + edge_marker.colors[i * 2].a = 1.0; + edge_marker.colors[i * 2 + 1].r = 1.0 - p2; + edge_marker.colors[i * 2 + 1].g = p2; + edge_marker.colors[i * 2 + 1].a = 1.0; + + + if(std::abs(v1->id() - v2->id()) > 2) { + edge_marker.points[i * 2].z += 0.5; + edge_marker.points[i * 2 + 1].z += 0.5; + } + + continue; + } + + g2o::EdgeSE3Plane* edge_plane = dynamic_cast(edge); + if(edge_plane) { + g2o::VertexSE3* v1 = dynamic_cast(edge_plane->vertices()[0]); + Eigen::Vector3d pt1 = v1->estimate().translation(); + Eigen::Vector3d pt2(pt1.x(), pt1.y(), 0.0); + + edge_marker.points[i * 2].x = pt1.x(); + edge_marker.points[i * 2].y = pt1.y(); + edge_marker.points[i * 2].z = pt1.z(); + edge_marker.points[i * 2 + 1].x = pt2.x(); + edge_marker.points[i * 2 + 1].y = pt2.y(); + edge_marker.points[i * 2 + 1].z = pt2.z(); + + edge_marker.colors[i * 2].b = 1.0; + edge_marker.colors[i * 2].a = 1.0; + edge_marker.colors[i * 2 + 1].b = 1.0; + edge_marker.colors[i * 2 + 1].a = 1.0; + + continue; + } + + g2o::EdgeSE3PriorXY* edge_priori_xy = dynamic_cast(edge); + if(edge_priori_xy) { + g2o::VertexSE3* v1 = dynamic_cast(edge_priori_xy->vertices()[0]); + Eigen::Vector3d pt1 = v1->estimate().translation(); + Eigen::Vector3d pt2 = Eigen::Vector3d::Zero(); + pt2.head<2>() = edge_priori_xy->measurement(); + + edge_marker.points[i * 2].x = pt1.x(); + edge_marker.points[i * 2].y = pt1.y(); + edge_marker.points[i * 2].z = pt1.z() + 0.5; + edge_marker.points[i * 2 + 1].x = pt2.x(); + edge_marker.points[i * 2 + 1].y = pt2.y(); + edge_marker.points[i * 2 + 1].z = pt2.z() + 0.5; + + edge_marker.colors[i * 2].r = 1.0; + edge_marker.colors[i * 2].a = 1.0; + edge_marker.colors[i * 2 + 1].r = 1.0; + edge_marker.colors[i * 2 + 1].a = 1.0; + + continue; + } + + g2o::EdgeSE3PriorXYZ* edge_priori_xyz = dynamic_cast(edge); + if(edge_priori_xyz) { + g2o::VertexSE3* v1 = dynamic_cast(edge_priori_xyz->vertices()[0]); + Eigen::Vector3d pt1 = v1->estimate().translation(); + Eigen::Vector3d pt2 = edge_priori_xyz->measurement(); + + edge_marker.points[i * 2].x = pt1.x(); + edge_marker.points[i * 2].y = pt1.y(); + edge_marker.points[i * 2].z = pt1.z() + 0.5; + edge_marker.points[i * 2 + 1].x = pt2.x(); + edge_marker.points[i * 2 + 1].y = pt2.y(); + edge_marker.points[i * 2 + 1].z = pt2.z(); + + edge_marker.colors[i * 2].r = 1.0; + edge_marker.colors[i * 2].a = 1.0; + edge_marker.colors[i * 2 + 1].r = 1.0; + edge_marker.colors[i * 2 + 1].a = 1.0; + + continue; + } + } + + // sphere + visualization_msgs::Marker& sphere_marker = markers.markers[3]; + sphere_marker.header.frame_id = "map"; + sphere_marker.header.stamp = stamp; + sphere_marker.ns = "loop_close_radius"; + sphere_marker.id = 3; + sphere_marker.type = visualization_msgs::Marker::SPHERE; + + if(!keyframes.empty()) { + Eigen::Vector3d pos = keyframes.back()->node->estimate().translation(); + sphere_marker.pose.position.x = pos.x(); + sphere_marker.pose.position.y = pos.y(); + sphere_marker.pose.position.z = pos.z(); + } + sphere_marker.pose.orientation.w = 1.0; + sphere_marker.scale.x = sphere_marker.scale.y = sphere_marker.scale.z = loop_detector->get_distance_thresh() * 2.0; + + sphere_marker.color.r = 1.0; + sphere_marker.color.a = 0.3; + + return markers; + } + + /** + * @brief dump all data to the current directory + * @param req + * @param res + * @return + */ + bool dump_service(hdl_graph_slam::DumpGraphRequest& req, hdl_graph_slam::DumpGraphResponse& res) { + std::lock_guard lock(main_thread_mutex); + + std::string directory = req.destination; + + if(directory.empty()) { + std::array buffer; + buffer.fill(0); + time_t rawtime; + time(&rawtime); + const auto timeinfo = localtime(&rawtime); + strftime(buffer.data(), sizeof(buffer), "%d-%m-%Y %H:%M:%S", timeinfo); + } + + if(!boost::filesystem::is_directory(directory)) { + boost::filesystem::create_directory(directory); + } + + std::cout << "all data dumped to:" << directory << std::endl; + + graph_slam->save(directory + "/graph.g2o"); + for(int i = 0; i < keyframes.size(); i++) { + std::stringstream sst; + sst << boost::format("%s/%06d") % directory % i; + + keyframes[i]->save(sst.str()); + } + + if(zero_utm) { + std::ofstream zero_utm_ofs(directory + "/zero_utm"); + zero_utm_ofs << boost::format("%.6f %.6f %.6f") % zero_utm->x() % zero_utm->y() % zero_utm->z() << std::endl; + } + + std::ofstream ofs(directory + "/special_nodes.csv"); + ofs << "anchor_node " << (anchor_node == nullptr ? -1 : anchor_node->id()) << std::endl; + ofs << "anchor_edge " << (anchor_edge == nullptr ? -1 : anchor_edge->id()) << std::endl; + ofs << "floor_node " << (floor_plane_node == nullptr ? -1 : floor_plane_node->id()) << std::endl; + + res.success = true; + return true; + } + + /** + * @brief save map data as pcd + * @param req + * @param res + * @return + */ + bool save_map_service(hdl_graph_slam::SaveMapRequest& req, hdl_graph_slam::SaveMapResponse& res) { + std::vector snapshot; + + keyframes_snapshot_mutex.lock(); + snapshot = keyframes_snapshot; + keyframes_snapshot_mutex.unlock(); + + auto cloud = map_cloud_generator->generate(snapshot, req.resolution); + if(!cloud) { + res.success = false; + return true; + } + + if(zero_utm && req.utm) { + for(auto& pt : cloud->points) { + pt.getVector3fMap() += (*zero_utm).cast(); + } + } + + cloud->header.frame_id = map_frame_id; + cloud->header.stamp = snapshot.back()->cloud->header.stamp; + + if(zero_utm) { + std::ofstream ofs(req.destination + ".utm"); + ofs << boost::format("%.6f %.6f %.6f") % zero_utm->x() % zero_utm->y() % zero_utm->z() << std::endl; + } + + int ret = pcl::io::savePCDFileBinary(req.destination, *cloud); + res.success = ret == 0; + + return true; + } + +private: + // ROS + ros::NodeHandle nh; + ros::NodeHandle mt_nh; + ros::NodeHandle private_nh; + ros::WallTimer optimization_timer; + ros::WallTimer map_publish_timer; + + std::unique_ptr> odom_sub; + std::unique_ptr> cloud_sub; + std::unique_ptr> sync; + + ros::Subscriber gps_sub; + ros::Subscriber nmea_sub; + ros::Subscriber navsat_sub; + + ros::Subscriber imu_sub; + ros::Subscriber floor_sub; + + ros::Publisher markers_pub; + + std::string map_frame_id; + std::string odom_frame_id; + + std::mutex trans_odom2map_mutex; + Eigen::Matrix4f trans_odom2map; + ros::Publisher odom2map_pub; + + std::string points_topic; + ros::Publisher read_until_pub; + ros::Publisher map_points_pub; + + tf::TransformListener tf_listener; + + ros::ServiceServer dump_service_server; + ros::ServiceServer save_map_service_server; + + // keyframe queue + std::string base_frame_id; + std::mutex keyframe_queue_mutex; + std::deque keyframe_queue; + + // gps queue + double gps_time_offset; + double gps_edge_stddev_xy; + double gps_edge_stddev_z; + boost::optional zero_utm; + std::mutex gps_queue_mutex; + std::deque gps_queue; + + // imu queue + double imu_time_offset; + bool enable_imu_orientation; + double imu_orientation_edge_stddev; + bool enable_imu_acceleration; + double imu_acceleration_edge_stddev; + std::mutex imu_queue_mutex; + std::deque imu_queue; + + // floor_coeffs queue + double floor_edge_stddev; + std::mutex floor_coeffs_queue_mutex; + std::deque floor_coeffs_queue; + + // for map cloud generation + std::atomic_bool graph_updated; + double map_cloud_resolution; + std::mutex keyframes_snapshot_mutex; + std::vector keyframes_snapshot; + std::unique_ptr map_cloud_generator; + + // graph slam + // all the below members must be accessed after locking main_thread_mutex + std::mutex main_thread_mutex; + + int max_keyframes_per_update; + std::deque new_keyframes; + + g2o::VertexSE3* anchor_node; + g2o::EdgeSE3* anchor_edge; + g2o::VertexPlane* floor_plane_node; + std::vector keyframes; + std::unordered_map keyframe_hash; + + std::unique_ptr graph_slam; + std::unique_ptr loop_detector; + std::unique_ptr keyframe_updater; + std::unique_ptr nmea_parser; + + std::string seq_; + bool use_test_; + + + std::unique_ptr inf_calclator; +}; + +} // namespace hdl_graph_slam + +PLUGINLIB_EXPORT_CLASS(hdl_graph_slam::HdlGraphSlamNodelet, nodelet::Nodelet) diff --git a/apps/prefiltering_nodelet.cpp b/apps/prefiltering_nodelet.cpp new file mode 100644 index 0000000..07bb648 --- /dev/null +++ b/apps/prefiltering_nodelet.cpp @@ -0,0 +1,270 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace hdl_graph_slam { + +class PrefilteringNodelet : public nodelet::Nodelet { +public: + typedef pcl::PointXYZI PointT; + + PrefilteringNodelet() {} + virtual ~PrefilteringNodelet() {} + + virtual void onInit() { + nh = getNodeHandle(); + private_nh = getPrivateNodeHandle(); + + initialize_params(); + + if(private_nh.param("deskewing", false)) { + imu_sub = nh.subscribe("/imu/data", 1, &PrefilteringNodelet::imu_callback, this); + } + + points_sub = nh.subscribe("/velodyne_points", 64, &PrefilteringNodelet::cloud_callback, this); + points_pub = nh.advertise("/filtered_points", 32); + colored_pub = nh.advertise("/colored_points", 32); + } + +private: + void initialize_params() { + std::string downsample_method = private_nh.param("downsample_method", "VOXELGRID"); + double downsample_resolution = private_nh.param("downsample_resolution", 0.1); + + if(downsample_method == "VOXELGRID") { + std::cout << "downsample: VOXELGRID " << downsample_resolution << std::endl; + auto voxelgrid = new pcl::VoxelGrid(); + voxelgrid->setLeafSize(downsample_resolution, downsample_resolution, downsample_resolution); + downsample_filter.reset(voxelgrid); + } else if(downsample_method == "APPROX_VOXELGRID") { + std::cout << "downsample: APPROX_VOXELGRID " << downsample_resolution << std::endl; + pcl::ApproximateVoxelGrid::Ptr approx_voxelgrid(new pcl::ApproximateVoxelGrid()); + approx_voxelgrid->setLeafSize(downsample_resolution, downsample_resolution, downsample_resolution); + downsample_filter = approx_voxelgrid; + } else { + if(downsample_method != "NONE") { + std::cerr << "warning: unknown downsampling type (" << downsample_method << ")" << std::endl; + std::cerr << " : use passthrough filter" << std::endl; + } + std::cout << "downsample: NONE" << std::endl; + } + + std::string outlier_removal_method = private_nh.param("outlier_removal_method", "STATISTICAL"); + if(outlier_removal_method == "STATISTICAL") { + int mean_k = private_nh.param("statistical_mean_k", 20); + double stddev_mul_thresh = private_nh.param("statistical_stddev", 1.0); + std::cout << "outlier_removal: STATISTICAL " << mean_k << " - " << stddev_mul_thresh << std::endl; + + pcl::StatisticalOutlierRemoval::Ptr sor(new pcl::StatisticalOutlierRemoval()); + sor->setMeanK(mean_k); + sor->setStddevMulThresh(stddev_mul_thresh); + outlier_removal_filter = sor; + } else if(outlier_removal_method == "RADIUS") { + double radius = private_nh.param("radius_radius", 0.8); + int min_neighbors = private_nh.param("radius_min_neighbors", 2); + std::cout << "outlier_removal: RADIUS " << radius << " - " << min_neighbors << std::endl; + + pcl::RadiusOutlierRemoval::Ptr rad(new pcl::RadiusOutlierRemoval()); + rad->setRadiusSearch(radius); + rad->setMinNeighborsInRadius(min_neighbors); + outlier_removal_filter = rad; + } else { + std::cout << "outlier_removal: NONE" << std::endl; + } + + use_distance_filter = private_nh.param("use_distance_filter", true); + distance_near_thresh = private_nh.param("distance_near_thresh", 1.0); + distance_far_thresh = private_nh.param("distance_far_thresh", 100.0); + + base_link_frame = private_nh.param("base_link_frame", ""); + } + + void imu_callback(const sensor_msgs::ImuConstPtr& imu_msg) { + imu_queue.push_back(imu_msg); + } + + void cloud_callback(const pcl::PointCloud& src_cloud_r) { + pcl::PointCloud::ConstPtr src_cloud = src_cloud_r.makeShared(); + if(src_cloud->empty()) { + return; + } + src_cloud = deskewing(src_cloud); + + // if base_link_frame is defined, transform the input cloud to the frame + if(!base_link_frame.empty()) { + if(!tf_listener.canTransform(base_link_frame, src_cloud->header.frame_id, ros::Time(0))) { + std::cerr << "failed to find transform between " << base_link_frame << " and " << src_cloud->header.frame_id << std::endl; + } + + tf::StampedTransform transform; + tf_listener.waitForTransform(base_link_frame, src_cloud->header.frame_id, ros::Time(0), ros::Duration(2.0)); + tf_listener.lookupTransform(base_link_frame, src_cloud->header.frame_id, ros::Time(0), transform); + + pcl::PointCloud::Ptr transformed(new pcl::PointCloud()); + pcl_ros::transformPointCloud(*src_cloud, *transformed, transform); + transformed->header.frame_id = base_link_frame; + transformed->header.stamp = src_cloud->header.stamp; + src_cloud = transformed; + } + + pcl::PointCloud::ConstPtr filtered = distance_filter(src_cloud); + filtered = downsample(filtered); + filtered = outlier_removal(filtered); + //std::cout<<"filtered cloud size: "<size()<::ConstPtr downsample(const pcl::PointCloud::ConstPtr& cloud) const { + if(!downsample_filter) { + return cloud; + } + + pcl::PointCloud::Ptr filtered(new pcl::PointCloud()); + downsample_filter->setInputCloud(cloud); + downsample_filter->filter(*filtered); + filtered->header = cloud->header; + + return filtered; + } + + pcl::PointCloud::ConstPtr outlier_removal(const pcl::PointCloud::ConstPtr& cloud) const { + if(!outlier_removal_filter) { + return cloud; + } + + pcl::PointCloud::Ptr filtered(new pcl::PointCloud()); + outlier_removal_filter->setInputCloud(cloud); + outlier_removal_filter->filter(*filtered); + filtered->header = cloud->header; + + return filtered; + } + + pcl::PointCloud::ConstPtr distance_filter(const pcl::PointCloud::ConstPtr& cloud) const { + pcl::PointCloud::Ptr filtered(new pcl::PointCloud()); + filtered->reserve(cloud->size()); + + std::copy_if(cloud->begin(), cloud->end(), std::back_inserter(filtered->points), [&](const PointT& p) { + double d = p.getVector3fMap().norm(); + return d > distance_near_thresh && d < distance_far_thresh; + }); + + filtered->width = filtered->size(); + filtered->height = 1; + filtered->is_dense = false; + + filtered->header = cloud->header; + + return filtered; + } + + pcl::PointCloud::ConstPtr deskewing(const pcl::PointCloud::ConstPtr& cloud) { + ros::Time stamp = pcl_conversions::fromPCL(cloud->header.stamp); + if(imu_queue.empty()) { + return cloud; + } + + // the color encodes the point number in the point sequence + if(colored_pub.getNumSubscribers()) { + pcl::PointCloud::Ptr colored(new pcl::PointCloud()); + colored->header = cloud->header; + colored->is_dense = cloud->is_dense; + colored->width = cloud->width; + colored->height = cloud->height; + colored->resize(cloud->size()); + + for(int i = 0; i < cloud->size(); i++) { + double t = static_cast(i) / cloud->size(); + colored->at(i).getVector4fMap() = cloud->at(i).getVector4fMap(); + colored->at(i).r = 255 * t; + colored->at(i).g = 128; + colored->at(i).b = 255 * (1 - t); + } + colored_pub.publish(*colored); + } + + sensor_msgs::ImuConstPtr imu_msg = imu_queue.front(); + + auto loc = imu_queue.begin(); + for(; loc != imu_queue.end(); loc++) { + imu_msg = (*loc); + if((*loc)->header.stamp > stamp) { + break; + } + } + + imu_queue.erase(imu_queue.begin(), loc); + + Eigen::Vector3f ang_v(imu_msg->angular_velocity.x, imu_msg->angular_velocity.y, imu_msg->angular_velocity.z); + ang_v *= -1; + + pcl::PointCloud::Ptr deskewed(new pcl::PointCloud()); + deskewed->header = cloud->header; + deskewed->is_dense = cloud->is_dense; + deskewed->width = cloud->width; + deskewed->height = cloud->height; + deskewed->resize(cloud->size()); + + double scan_period = private_nh.param("scan_period", 0.1); + for(int i = 0; i < cloud->size(); i++) { + const auto& pt = cloud->at(i); + + // TODO: transform IMU data into the LIDAR frame + double delta_t = scan_period * static_cast(i) / cloud->size(); + Eigen::Quaternionf delta_q(1, delta_t / 2.0 * ang_v[0], delta_t / 2.0 * ang_v[1], delta_t / 2.0 * ang_v[2]); + Eigen::Vector3f pt_ = delta_q.inverse() * pt.getVector3fMap(); + + deskewed->at(i) = cloud->at(i); + deskewed->at(i).getVector3fMap() = pt_; + } + + return deskewed; + } + +private: + ros::NodeHandle nh; + ros::NodeHandle private_nh; + + ros::Subscriber imu_sub; + std::vector imu_queue; + + ros::Subscriber points_sub; + ros::Publisher points_pub; + + ros::Publisher colored_pub; + + tf::TransformListener tf_listener; + + std::string base_link_frame; + + bool use_distance_filter; + double distance_near_thresh; + double distance_far_thresh; + + pcl::Filter::Ptr downsample_filter; + pcl::Filter::Ptr outlier_removal_filter; +}; + +} // namespace hdl_graph_slam + +PLUGINLIB_EXPORT_CLASS(hdl_graph_slam::PrefilteringNodelet, nodelet::Nodelet) diff --git a/apps/scan_matching_odometry_nodelet.cpp b/apps/scan_matching_odometry_nodelet.cpp new file mode 100644 index 0000000..9480d65 --- /dev/null +++ b/apps/scan_matching_odometry_nodelet.cpp @@ -0,0 +1,430 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace hdl_graph_slam { + +class ScanMatchingOdometryNodelet : public nodelet::Nodelet { +public: + FeatureAssociation FA; + typedef pcl::PointXYZI PointT; + EIGEN_MAKE_ALIGNED_OPERATOR_NEW; + + ScanMatchingOdometryNodelet() {} + virtual ~ScanMatchingOdometryNodelet() {} + + virtual void onInit() { + NODELET_DEBUG("initializing scan_matching_odometry_nodelet..."); + nh = getNodeHandle(); + private_nh = getPrivateNodeHandle(); + + initialize_params(); + + if(private_nh.param("enable_imu_frontend", false)) { + msf_pose_sub = nh.subscribe("/msf_core/pose", 1, boost::bind(&ScanMatchingOdometryNodelet::msf_pose_callback, this, _1, false)); + msf_pose_after_update_sub = nh.subscribe("/msf_core/pose_after_update", 1, boost::bind(&ScanMatchingOdometryNodelet::msf_pose_callback, this, _1, true)); + } + + points_sub = nh.subscribe("/filtered_points", 256, &ScanMatchingOdometryNodelet::cloud_callback, this); + read_until_pub = nh.advertise("/scan_matching_odometry/read_until", 32); + odom_pub = nh.advertise("/odom", 32); + trans_pub = nh.advertise("/scan_matching_odometry/transform", 32); + status_pub = private_nh.advertise("/scan_matching_odometry/status", 8); + aligned_points_pub = nh.advertise("/aligned_points", 32); + } + +private: + /** + * @brief initialize parameters + */ + void initialize_params() { + auto& pnh = private_nh; + points_topic = pnh.param("points_topic", "/velodyne_points"); + odom_frame_id = pnh.param("odom_frame_id", "odom"); + robot_odom_frame_id = pnh.param("robot_odom_frame_id", "robot_odom"); + + // The minimum tranlational distance and rotation angle between keyframes. + // If this value is zero, frames are always compared with the previous frame + keyframe_delta_trans = pnh.param("keyframe_delta_trans", 0.25); + keyframe_delta_angle = pnh.param("keyframe_delta_angle", 0.15); + keyframe_delta_time = pnh.param("keyframe_delta_time", 1.0); + + Thre_ = pnh.param("Thre_", 1.0); + Time_ = pnh.param("Time_", 1.0); + std::cout<<"Thre_ and Time_: "<("transform_thresholding", false); + max_acceptable_trans = pnh.param("max_acceptable_trans", 1.0); + max_acceptable_angle = pnh.param("max_acceptable_angle", 1.0); + + // select a downsample method (VOXELGRID, APPROX_VOXELGRID, NONE) + std::string downsample_method = pnh.param("downsample_method", "VOXELGRID"); + double downsample_resolution = pnh.param("downsample_resolution", 0.25); + if(downsample_method == "VOXELGRID") { + std::cout << "downsample: VOXELGRID " << downsample_resolution << std::endl; + auto voxelgrid = new pcl::VoxelGrid(); + voxelgrid->setLeafSize(downsample_resolution, downsample_resolution, downsample_resolution); + downsample_filter.reset(voxelgrid); + } else if(downsample_method == "APPROX_VOXELGRID") { + std::cout << "downsample: APPROX_VOXELGRID " << downsample_resolution << std::endl; + pcl::ApproximateVoxelGrid::Ptr approx_voxelgrid(new pcl::ApproximateVoxelGrid()); + approx_voxelgrid->setLeafSize(downsample_resolution, downsample_resolution, downsample_resolution); + downsample_filter = approx_voxelgrid; + } else { + if(downsample_method != "NONE") { + std::cerr << "warning: unknown downsampling type (" << downsample_method << ")" << std::endl; + std::cerr << " : use passthrough filter" << std::endl; + } + std::cout << "downsample: NONE" << std::endl; + pcl::PassThrough::Ptr passthrough(new pcl::PassThrough()); + downsample_filter = passthrough; + } + registration = select_registration_method(pnh); + + + } + + /** + * @brief callback for point clouds + * @param cloud_msg point cloud msg + */ + void cloud_callback(const sensor_msgs::PointCloud2ConstPtr& cloud_msg) { + if(!ros::ok()) { + return; + } + + pcl::PointCloud::Ptr cloud(new pcl::PointCloud()); + pcl::fromROSMsg(*cloud_msg, *cloud); + + Eigen::Matrix4f pose = matching(cloud_msg->header.stamp, cloud); + publish_odometry(cloud_msg->header.stamp, cloud_msg->header.frame_id, pose); + + // In offline estimation, point clouds until the published time will be supplied + std_msgs::HeaderPtr read_until(new std_msgs::Header()); + read_until->frame_id = points_topic; + read_until->stamp = cloud_msg->header.stamp + ros::Duration(1, 0); + read_until_pub.publish(read_until); + + read_until->frame_id = "/filtered_points"; + read_until_pub.publish(read_until); + } + + void msf_pose_callback(const geometry_msgs::PoseWithCovarianceStampedConstPtr& pose_msg, bool after_update) { + if(after_update) { + msf_pose_after_update = pose_msg; + } else { + msf_pose = pose_msg; + } + } + + /** + * @brief downsample a point cloud + * @param cloud input cloud + * @return downsampled point cloud + */ + pcl::PointCloud::ConstPtr downsample(const pcl::PointCloud::ConstPtr& cloud) const { + if(!downsample_filter) { + return cloud; + } + + pcl::PointCloud::Ptr filtered(new pcl::PointCloud()); + downsample_filter->setInputCloud(cloud); + downsample_filter->filter(*filtered); + + return filtered; + } + + /** + * @brief estimate the relative pose between an input cloud and a keyframe cloud + * @param stamp the timestamp of the input cloud + * @param cloud the input cloud + * @return the relative pose between the input cloud and the keyframe cloud + */ + Eigen::Matrix4f matching(const ros::Time& stamp, const pcl::PointCloud::ConstPtr& cloud) { + if(!keyframe) { + prev_time = ros::Time(); + prev_trans.setIdentity(); + prev_trans0_.setIdentity(); + keyframe_pose.setIdentity(); + keyframe_stamp = stamp; +// + keyframe = downsample(cloud); + registration->setInputTarget(keyframe); + Accumulated_trans_.setIdentity(); + FA.setInitTargetCloud(cloud); + last_frame_ = cloud; + return Eigen::Matrix4f::Identity(); + } + frames_num_++; + Eigen::Matrix4f trans0; + clock_t start_s = clock(); + FA.setSourceCloud(cloud); + clock_t end_s = clock(); + std::cout<<"setSourceCloud "<<(double)(end_s-start_s)/CLOCKS_PER_SEC<("enable_imu_frontend", false)) { + if(msf_pose && msf_pose->header.stamp > keyframe_stamp && msf_pose_after_update && msf_pose_after_update->header.stamp > keyframe_stamp) { + Eigen::Isometry3d pose0 = pose2isometry(msf_pose_after_update->pose.pose); + Eigen::Isometry3d pose1 = pose2isometry(msf_pose->pose.pose); + Eigen::Isometry3d delta = pose0.inverse() * pose1; + msf_source = "imu"; + msf_delta = delta.cast(); + } else { + std::cerr << "msf data is too old" << std::endl; + } + } else if(private_nh.param("enable_robot_odometry_init_guess", false) && !prev_time.isZero()) { + tf::StampedTransform transform; + if(tf_listener.waitForTransform(cloud->header.frame_id, stamp, cloud->header.frame_id, prev_time, robot_odom_frame_id, ros::Duration(0))) { + tf_listener.lookupTransform(cloud->header.frame_id, stamp, cloud->header.frame_id, prev_time, robot_odom_frame_id, transform); + } else if(tf_listener.waitForTransform(cloud->header.frame_id, ros::Time(0), cloud->header.frame_id, prev_time, robot_odom_frame_id, ros::Duration(0))) { + tf_listener.lookupTransform(cloud->header.frame_id, ros::Time(0), cloud->header.frame_id, prev_time, robot_odom_frame_id, transform); + } + + if(transform.stamp_.isZero()) { + NODELET_WARN_STREAM("failed to look up transform between " << cloud->header.frame_id << " and " << robot_odom_frame_id); + } else { + msf_source = "odometry"; + msf_delta = tf2isometry(transform).cast(); + } + } + + //publish_scan_matching_status(stamp, cloud->header.frame_id, aligned, msf_source, msf_delta); + trans0 = FA.getFinalTransformation(); + //std::cout<<"trans0: "<(0, 3).norm(); + double da = std::acos(Eigen::Quaternionf(delta.block<3, 3>(0, 0)).w()); + //if uniform motion is not satisfied, will not accept this result. + if(dx > max_acceptable_trans || da > max_acceptable_angle) { + NODELET_INFO_STREAM("too large transform!! " << dx << "[m] " << da << "[rad]"); + NODELET_INFO_STREAM("ignore this frame(" << stamp << ")"); + return keyframe_pose * prev_trans; + } + } + prev_trans0_ = trans0; + prev_time = stamp; + prev_trans = trans; + + auto keyframe_trans = matrix2transform(stamp, keyframe_pose, odom_frame_id, "keyframe"); + keyframe_broadcaster.sendTransform(keyframe_trans); + + double delta_trans = trans.block<3, 1>(0, 3).norm(); + double delta_angle = std::acos(Eigen::Quaternionf(trans.block<3, 3>(0, 0)).w()); + double delta_time = (stamp - keyframe_stamp).toSec(); + pcl::PointCloud::Ptr aligned(new pcl::PointCloud()); + // + std::cout<<"delta_trans and time: "< keyframe_delta_trans || delta_angle > keyframe_delta_angle || delta_time > keyframe_delta_time) { + registration->setInputSource(filtered); + clock_t start = clock(); + registration->align(*aligned, prev_trans * msf_delta.matrix()); + clock_t end = clock(); + std::cout<<"GICP running speed "<<(double)(end-start)/CLOCKS_PER_SEC<hasConverged()) { + NODELET_INFO_STREAM("scan matching has not converged!!"); + NODELET_INFO_STREAM("ignore this frame(" << stamp << ")"); + return keyframe_pose * prev_trans; + } + keyframe = filtered; + registration->setInputTarget(keyframe); + Accumulated_trans_.setIdentity(); + odom = keyframe_pose * registration->getFinalTransformation(); + keyframe_pose = odom; + keyframe_stamp = stamp; + prev_time = stamp; + prev_trans.setIdentity(); + }} + + if (aligned_points_pub.getNumSubscribers() > 0) + { + pcl::transformPointCloud (*cloud, *aligned, odom); + aligned->header.frame_id=odom_frame_id; + aligned_points_pub.publish(*aligned); + } + return odom; + } + + /** + * @brief publish odometry + * @param stamp timestamp + * @param pose odometry pose to be published + */ + void publish_odometry(const ros::Time& stamp, const std::string& base_frame_id, const Eigen::Matrix4f& pose) { + // publish transform stamped for IMU integration + geometry_msgs::TransformStamped odom_trans = matrix2transform(stamp, pose, odom_frame_id, base_frame_id); + trans_pub.publish(odom_trans); + + // broadcast the transform over tf + odom_broadcaster.sendTransform(odom_trans); + + // publish the transform + nav_msgs::Odometry odom; + odom.header.stamp = stamp; + odom.header.frame_id = odom_frame_id; + + odom.pose.pose.position.x = pose(0, 3); + odom.pose.pose.position.y = pose(1, 3); + odom.pose.pose.position.z = pose(2, 3); + odom.pose.pose.orientation = odom_trans.transform.rotation; + + odom.child_frame_id = base_frame_id; + odom.twist.twist.linear.x = 0.0; + odom.twist.twist.linear.y = 0.0; + odom.twist.twist.angular.z = 0.0; + + odom_pub.publish(odom); + } + + /** + * @brief publish scan matching status + */ + void publish_scan_matching_status(const ros::Time& stamp, const std::string& frame_id, pcl::PointCloud::ConstPtr aligned, const std::string& msf_source, const Eigen::Isometry3f& msf_delta) { + if(!status_pub.getNumSubscribers()) { + return; + } + + ScanMatchingStatus status; + status.header.frame_id = frame_id; + status.header.stamp = stamp; + status.has_converged = registration->hasConverged(); + status.matching_error = registration->getFitnessScore(); + + const double max_correspondence_dist = 0.5; + + int num_inliers = 0; + std::vector k_indices; + std::vector k_sq_dists; + for(int i=0; isize(); i++) { + const auto& pt = aligned->at(i); + registration->getSearchMethodTarget()->nearestKSearch(pt, 1, k_indices, k_sq_dists); + if(k_sq_dists[0] < max_correspondence_dist * max_correspondence_dist) { + num_inliers++; + } + } + status.inlier_fraction = static_cast(num_inliers) / aligned->size(); + + status.relative_pose = isometry2pose(Eigen::Isometry3f(registration->getFinalTransformation()).cast()); + + if(!msf_source.empty()) { + status.prediction_labels.resize(1); + status.prediction_labels[0].data = msf_source; + + status.prediction_errors.resize(1); + Eigen::Isometry3f error = Eigen::Isometry3f(registration->getFinalTransformation()).inverse() * msf_delta; + status.prediction_errors[0] = isometry2pose(error.cast()); + } + + status_pub.publish(status); + } + +private: + // ROS topics + ros::NodeHandle nh; + ros::NodeHandle private_nh; + + ros::Subscriber points_sub; + ros::Subscriber msf_pose_sub; + ros::Subscriber msf_pose_after_update_sub; + + ros::Publisher odom_pub; + ros::Publisher trans_pub; + ros::Publisher aligned_points_pub; + ros::Publisher status_pub; + tf::TransformListener tf_listener; + tf::TransformBroadcaster odom_broadcaster; + tf::TransformBroadcaster keyframe_broadcaster; + + std::string points_topic; + std::string odom_frame_id; + std::string robot_odom_frame_id; + ros::Publisher read_until_pub; + + // keyframe parameters + double keyframe_delta_trans; // minimum distance between keyframes + double keyframe_delta_angle; // + double keyframe_delta_time; // + + // registration validation by thresholding + bool transform_thresholding; // + double max_acceptable_trans; // + double max_acceptable_angle; + + // odometry calculation + geometry_msgs::PoseWithCovarianceStampedConstPtr msf_pose; + geometry_msgs::PoseWithCovarianceStampedConstPtr msf_pose_after_update; + + ros::Time prev_time; + Eigen::Matrix4f prev_trans; // previous estimated transform from keyframe + Eigen::Matrix4f keyframe_pose; // keyframe pose + ros::Time keyframe_stamp; // keyframe time + pcl::PointCloud::ConstPtr keyframe; // keyframe point cloud + pcl::PointCloud::ConstPtr last_frame_; // keyframe point cloud + Eigen::Matrix4f Accumulated_trans_; + Eigen::Matrix4f prev_trans0_; //used for save the scan-to-scan matching + pcl::PointCloud Accumulated_frames_; + int Accumulated_no_; + int frames_num_; + + double Time_; + double Thre_; + + // + pcl::Filter::Ptr downsample_filter; + pcl::Registration::Ptr registration; +}; + +} // namespace hdl_graph_slam + +PLUGINLIB_EXPORT_CLASS(hdl_graph_slam::ScanMatchingOdometryNodelet, nodelet::Nodelet) diff --git a/apps/scan_matching_odometry_nodelet_raw.cpp b/apps/scan_matching_odometry_nodelet_raw.cpp new file mode 100644 index 0000000..fb95159 --- /dev/null +++ b/apps/scan_matching_odometry_nodelet_raw.cpp @@ -0,0 +1,399 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include + +namespace hdl_graph_slam { + +class ScanMatchingOdometryNodelet : public nodelet::Nodelet { +public: + typedef pcl::PointXYZI PointT; + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + + ScanMatchingOdometryNodelet() {} + virtual ~ScanMatchingOdometryNodelet() {} + + virtual void onInit() { + NODELET_DEBUG("initializing scan_matching_odometry_nodelet..."); + nh = getNodeHandle(); + private_nh = getPrivateNodeHandle(); + + initialize_params(); + + if(private_nh.param("enable_imu_frontend", false)) { + msf_pose_sub = nh.subscribe("/msf_core/pose", 1, boost::bind(&ScanMatchingOdometryNodelet::msf_pose_callback, this, _1, false)); + msf_pose_after_update_sub = nh.subscribe("/msf_core/pose_after_update", 1, boost::bind(&ScanMatchingOdometryNodelet::msf_pose_callback, this, _1, true)); + } + + points_sub = nh.subscribe("/filtered_points", 256, &ScanMatchingOdometryNodelet::cloud_callback, this); + read_until_pub = nh.advertise("/scan_matching_odometry/read_until", 32); + odom_pub = nh.advertise("/odom", 32); + trans_pub = nh.advertise("/scan_matching_odometry/transform", 32); + status_pub = private_nh.advertise("/scan_matching_odometry/status", 8); + aligned_points_pub = nh.advertise("/aligned_points", 32); + } + +private: + /** + * @brief initialize parameters + */ + void initialize_params() { + auto& pnh = private_nh; + points_topic = pnh.param("points_topic", "/velodyne_points"); + odom_frame_id = pnh.param("odom_frame_id", "odom"); + robot_odom_frame_id = pnh.param("robot_odom_frame_id", "robot_odom"); + + // The minimum tranlational distance and rotation angle between keyframes. + // If this value is zero, frames are always compared with the previous frame + keyframe_delta_trans = pnh.param("keyframe_delta_trans", 0.25); + keyframe_delta_angle = pnh.param("keyframe_delta_angle", 0.15); + keyframe_delta_time = pnh.param("keyframe_delta_time", 1.0); + + // Registration validation by thresholding + transform_thresholding = pnh.param("transform_thresholding", false); + max_acceptable_trans = pnh.param("max_acceptable_trans", 1.0); + max_acceptable_angle = pnh.param("max_acceptable_angle", 1.0); + + // select a downsample method (VOXELGRID, APPROX_VOXELGRID, NONE) + std::string downsample_method = pnh.param("downsample_method", "VOXELGRID"); + double downsample_resolution = pnh.param("downsample_resolution", 0.1); + if(downsample_method == "VOXELGRID") { + std::cout << "downsample: VOXELGRID " << downsample_resolution << std::endl; + auto voxelgrid = new pcl::VoxelGrid(); + voxelgrid->setLeafSize(downsample_resolution, downsample_resolution, downsample_resolution); + downsample_filter.reset(voxelgrid); + } else if(downsample_method == "APPROX_VOXELGRID") { + std::cout << "downsample: APPROX_VOXELGRID " << downsample_resolution << std::endl; + pcl::ApproximateVoxelGrid::Ptr approx_voxelgrid(new pcl::ApproximateVoxelGrid()); + approx_voxelgrid->setLeafSize(downsample_resolution, downsample_resolution, downsample_resolution); + downsample_filter = approx_voxelgrid; + } else { + if(downsample_method != "NONE") { + std::cerr << "warning: unknown downsampling type (" << downsample_method << ")" << std::endl; + std::cerr << " : use passthrough filter" << std::endl; + } + std::cout << "downsample: NONE" << std::endl; + pcl::PassThrough::Ptr passthrough(new pcl::PassThrough()); + downsample_filter = passthrough; + } + + registration = select_registration_method(pnh); + } + + /** + * @brief callback for point clouds + * @param cloud_msg point cloud msg + */ + void cloud_callback(const sensor_msgs::PointCloud2ConstPtr& cloud_msg) { + if(!ros::ok()) { + return; + } + + pcl::PointCloud::Ptr cloud(new pcl::PointCloud()); + pcl::fromROSMsg(*cloud_msg, *cloud); + clock_t start = clock(); + Eigen::Matrix4f pose = matching(cloud_msg->header.stamp, cloud); + clock_t end = clock(); + std::cout<<"matching time cost: "<<(double)(end-start)/CLOCKS_PER_SEC<header.stamp, cloud_msg->header.frame_id, pose); + + // In offline estimation, point clouds until the published time will be supplied + std_msgs::HeaderPtr read_until(new std_msgs::Header()); + read_until->frame_id = points_topic; + read_until->stamp = cloud_msg->header.stamp + ros::Duration(1, 0); + read_until_pub.publish(read_until); + + read_until->frame_id = "/filtered_points"; + read_until_pub.publish(read_until); + } + + void msf_pose_callback(const geometry_msgs::PoseWithCovarianceStampedConstPtr& pose_msg, bool after_update) { + if(after_update) { + msf_pose_after_update = pose_msg; + } else { + msf_pose = pose_msg; + } + } + + /** + * @brief downsample a point cloud + * @param cloud input cloud + * @return downsampled point cloud + */ + pcl::PointCloud::ConstPtr downsample(const pcl::PointCloud::ConstPtr& cloud) const { + if(!downsample_filter) { + return cloud; + } + + pcl::PointCloud::Ptr filtered(new pcl::PointCloud()); + downsample_filter->setInputCloud(cloud); + downsample_filter->filter(*filtered); + + return filtered; + } + + /** + * @brief estimate the relative pose between an input cloud and a keyframe cloud + * @param stamp the timestamp of the input cloud + * @param cloud the input cloud + * @return the relative pose between the input cloud and the keyframe cloud + */ + Eigen::Matrix4f matching(const ros::Time& stamp, const pcl::PointCloud::ConstPtr& cloud) { + if(!keyframe) { + prev_time = ros::Time(); + prev_trans.setIdentity(); + keyframe_pose.setIdentity(); + keyframe_stamp = stamp; + keyframe = downsample(cloud); + registration->setInputTarget(keyframe); + frame_num_ = 0; + return Eigen::Matrix4f::Identity(); + } +// + auto filtered = downsample(cloud); + registration->setInputSource(filtered); + + std::string msf_source; + Eigen::Isometry3f msf_delta = Eigen::Isometry3f::Identity(); + + if(private_nh.param("enable_imu_frontend", false)) { + if(msf_pose && msf_pose->header.stamp > keyframe_stamp && msf_pose_after_update && msf_pose_after_update->header.stamp > keyframe_stamp) { + Eigen::Isometry3d pose0 = pose2isometry(msf_pose_after_update->pose.pose); + Eigen::Isometry3d pose1 = pose2isometry(msf_pose->pose.pose); + Eigen::Isometry3d delta = pose0.inverse() * pose1; + + msf_source = "imu"; + msf_delta = delta.cast(); + } else { + std::cerr << "msf data is too old" << std::endl; + } + } else if(private_nh.param("enable_robot_odometry_init_guess", false) && !prev_time.isZero()) { + tf::StampedTransform transform; + if(tf_listener.waitForTransform(cloud->header.frame_id, stamp, cloud->header.frame_id, prev_time, robot_odom_frame_id, ros::Duration(0))) { + tf_listener.lookupTransform(cloud->header.frame_id, stamp, cloud->header.frame_id, prev_time, robot_odom_frame_id, transform); + } else if(tf_listener.waitForTransform(cloud->header.frame_id, ros::Time(0), cloud->header.frame_id, prev_time, robot_odom_frame_id, ros::Duration(0))) { + tf_listener.lookupTransform(cloud->header.frame_id, ros::Time(0), cloud->header.frame_id, prev_time, robot_odom_frame_id, transform); + } + + if(transform.stamp_.isZero()) { + NODELET_WARN_STREAM("failed to look up transform between " << cloud->header.frame_id << " and " << robot_odom_frame_id); + } else { + msf_source = "odometry"; + msf_delta = tf2isometry(transform).cast(); + } + } + + pcl::PointCloud::Ptr aligned(new pcl::PointCloud()); + clock_t t0 = clock(); + registration->align(*aligned, prev_trans * msf_delta.matrix()); + clock_t t1 = clock(); + std::cout<<"align time: "<<(double)(t1-t0)/CLOCKS_PER_SEC<header.frame_id, aligned, msf_source, msf_delta); + clock_t t2 = clock(); + std::cout<<"publish_scan_matching_status time: "<<(double)(t2-t1)/CLOCKS_PER_SEC<hasConverged()) { + NODELET_INFO_STREAM("scan matching has not converged!!"); + NODELET_INFO_STREAM("ignore this frame(" << stamp << ")"); + return keyframe_pose * prev_trans; + } + + Eigen::Matrix4f trans = registration->getFinalTransformation(); + + //std::cout<<"trans: "<(0, 3).norm(); + double da = std::acos(Eigen::Quaternionf(delta.block<3, 3>(0, 0)).w()); + + if(dx > max_acceptable_trans || da > max_acceptable_angle) { + NODELET_INFO_STREAM("too large transform!! " << dx << "[m] " << da << "[rad]"); + NODELET_INFO_STREAM("ignore this frame(" << stamp << ")"); + return keyframe_pose * prev_trans; + } + } + + prev_time = stamp; + prev_trans = trans; + + auto keyframe_trans = matrix2transform(stamp, keyframe_pose, odom_frame_id, "keyframe"); + keyframe_broadcaster.sendTransform(keyframe_trans); + + double delta_trans = trans.block<3, 1>(0, 3).norm(); + double delta_angle = std::acos(Eigen::Quaternionf(trans.block<3, 3>(0, 0)).w()); + double delta_time = (stamp - keyframe_stamp).toSec(); + //std::cout<<"delta_time: "< keyframe_delta_trans || delta_angle > keyframe_delta_angle || delta_time > keyframe_delta_time) { + keyframe = filtered; + registration->setInputTarget(keyframe); + //use fastgicp + keyframe_pose = odom; + keyframe_stamp = stamp; + prev_time = stamp; + prev_trans.setIdentity(); + } + + if (aligned_points_pub.getNumSubscribers() > 0) + { + pcl::transformPointCloud(*cloud, *aligned, odom); + aligned->header.frame_id=odom_frame_id; + aligned_points_pub.publish(*aligned); + } + frame_num_++; + //if(frame_num_>100) + // sleep(1000); + + return odom; + } + + /** + * @brief publish odometry + * @param stamp timestamp + * @param pose odometry pose to be published + */ + void publish_odometry(const ros::Time& stamp, const std::string& base_frame_id, const Eigen::Matrix4f& pose) { + // publish transform stamped for IMU integration + geometry_msgs::TransformStamped odom_trans = matrix2transform(stamp, pose, odom_frame_id, base_frame_id); + trans_pub.publish(odom_trans); + + // broadcast the transform over tf + odom_broadcaster.sendTransform(odom_trans); + + // publish the transform + nav_msgs::Odometry odom; + odom.header.stamp = stamp; + odom.header.frame_id = odom_frame_id; + + odom.pose.pose.position.x = pose(0, 3); + odom.pose.pose.position.y = pose(1, 3); + odom.pose.pose.position.z = pose(2, 3); + odom.pose.pose.orientation = odom_trans.transform.rotation; + + odom.child_frame_id = base_frame_id; + odom.twist.twist.linear.x = 0.0; + odom.twist.twist.linear.y = 0.0; + odom.twist.twist.angular.z = 0.0; + + odom_pub.publish(odom); + } + + /** + * @brief publish scan matching status + */ + void publish_scan_matching_status(const ros::Time& stamp, const std::string& frame_id, pcl::PointCloud::ConstPtr aligned, const std::string& msf_source, const Eigen::Isometry3f& msf_delta) { + if(!status_pub.getNumSubscribers()) { + return; + } + + ScanMatchingStatus status; + status.header.frame_id = frame_id; + status.header.stamp = stamp; + status.has_converged = registration->hasConverged(); + status.matching_error = registration->getFitnessScore(); + + const double max_correspondence_dist = 0.5; + + int num_inliers = 0; + std::vector k_indices; + std::vector k_sq_dists; + for(int i=0; isize(); i++) { + const auto& pt = aligned->at(i); + registration->getSearchMethodTarget()->nearestKSearch(pt, 1, k_indices, k_sq_dists); + if(k_sq_dists[0] < max_correspondence_dist * max_correspondence_dist) { + num_inliers++; + } + } + status.inlier_fraction = static_cast(num_inliers) / aligned->size(); + + status.relative_pose = isometry2pose(Eigen::Isometry3f(registration->getFinalTransformation()).cast()); + + if(!msf_source.empty()) { + status.prediction_labels.resize(1); + status.prediction_labels[0].data = msf_source; + + status.prediction_errors.resize(1); + Eigen::Isometry3f error = Eigen::Isometry3f(registration->getFinalTransformation()).inverse() * msf_delta; + status.prediction_errors[0] = isometry2pose(error.cast()); + } + + status_pub.publish(status); + } + +private: + // ROS topics + ros::NodeHandle nh; + ros::NodeHandle private_nh; + + ros::Subscriber points_sub; + ros::Subscriber msf_pose_sub; + ros::Subscriber msf_pose_after_update_sub; + + ros::Publisher odom_pub; + ros::Publisher trans_pub; + ros::Publisher aligned_points_pub; + ros::Publisher status_pub; + tf::TransformListener tf_listener; + tf::TransformBroadcaster odom_broadcaster; + tf::TransformBroadcaster keyframe_broadcaster; + + std::string points_topic; + std::string odom_frame_id; + std::string robot_odom_frame_id; + ros::Publisher read_until_pub; + + // keyframe parameters + double keyframe_delta_trans; // minimum distance between keyframes + double keyframe_delta_angle; // + double keyframe_delta_time; // + + // registration validation by thresholding + bool transform_thresholding; // + double max_acceptable_trans; // + double max_acceptable_angle; + + // odometry calculation + geometry_msgs::PoseWithCovarianceStampedConstPtr msf_pose; + geometry_msgs::PoseWithCovarianceStampedConstPtr msf_pose_after_update; + + ros::Time prev_time; + Eigen::Matrix4f prev_trans; // previous estimated transform from keyframe + Eigen::Matrix4f keyframe_pose; // keyframe pose + ros::Time keyframe_stamp; // keyframe time + pcl::PointCloud::ConstPtr keyframe; // keyframe point cloud + int frame_num_; + + // + pcl::Filter::Ptr downsample_filter; + pcl::Registration::Ptr registration; +}; + +} // namespace hdl_graph_slam + +PLUGINLIB_EXPORT_CLASS(hdl_graph_slam::ScanMatchingOdometryNodelet, nodelet::Nodelet) diff --git a/cmake/FindG2O.cmake b/cmake/FindG2O.cmake new file mode 100644 index 0000000..a2438ca --- /dev/null +++ b/cmake/FindG2O.cmake @@ -0,0 +1,117 @@ +# Find the header files + +FIND_PATH(G2O_INCLUDE_DIR g2o/core/base_vertex.h + ${G2O_ROOT}/include + $ENV{G2O_ROOT}/include + $ENV{G2O_ROOT} + /usr/local/include + /usr/include + /opt/local/include + /sw/local/include + /sw/include + /opt/ros/$ENV{ROS_DISTRO}/include + NO_DEFAULT_PATH + ) + +# Macro to unify finding both the debug and release versions of the +# libraries; this is adapted from the OpenSceneGraph FIND_LIBRARY +# macro. + +MACRO(FIND_G2O_LIBRARY MYLIBRARY MYLIBRARYNAME) + + FIND_LIBRARY("${MYLIBRARY}_DEBUG" + NAMES "g2o_${MYLIBRARYNAME}_d" + PATHS + ${G2O_ROOT}/lib/Debug + ${G2O_ROOT}/lib + $ENV{G2O_ROOT}/lib/Debug + $ENV{G2O_ROOT}/lib + /opt/ros/$ENV{ROS_DISTRO}/lib + NO_DEFAULT_PATH + ) + + FIND_LIBRARY("${MYLIBRARY}_DEBUG" + NAMES "g2o_${MYLIBRARYNAME}_d" + PATHS + ~/Library/Frameworks + /Library/Frameworks + /usr/local/lib + /usr/local/lib64 + /usr/lib + /usr/lib64 + /opt/local/lib + /sw/local/lib + /sw/lib + ) + + FIND_LIBRARY(${MYLIBRARY} + NAMES "g2o_${MYLIBRARYNAME}" + PATHS + ${G2O_ROOT}/lib/Release + ${G2O_ROOT}/lib + $ENV{G2O_ROOT}/lib/Release + $ENV{G2O_ROOT}/lib + /opt/ros/$ENV{ROS_DISTRO}/lib + NO_DEFAULT_PATH + ) + + FIND_LIBRARY(${MYLIBRARY} + NAMES "g2o_${MYLIBRARYNAME}" + PATHS + ~/Library/Frameworks + /Library/Frameworks + /usr/local/lib + /usr/local/lib64 + /usr/lib + /usr/lib64 + /opt/local/lib + /sw/local/lib + /sw/lib + ) + + IF(NOT ${MYLIBRARY}_DEBUG) + IF(MYLIBRARY) + SET(${MYLIBRARY}_DEBUG ${MYLIBRARY}) + ENDIF(MYLIBRARY) + ENDIF( NOT ${MYLIBRARY}_DEBUG) + +ENDMACRO(FIND_G2O_LIBRARY LIBRARY LIBRARYNAME) + +# Find the core elements +FIND_G2O_LIBRARY(G2O_STUFF_LIBRARY stuff) +FIND_G2O_LIBRARY(G2O_CORE_LIBRARY core) + +# Find the CLI library +FIND_G2O_LIBRARY(G2O_CLI_LIBRARY cli) + +# Find the pluggable solvers +FIND_G2O_LIBRARY(G2O_SOLVER_CHOLMOD solver_cholmod) +FIND_G2O_LIBRARY(G2O_SOLVER_CSPARSE solver_csparse) +FIND_G2O_LIBRARY(G2O_SOLVER_CSPARSE_EXTENSION csparse_extension) +FIND_G2O_LIBRARY(G2O_SOLVER_DENSE solver_dense) +FIND_G2O_LIBRARY(G2O_SOLVER_PCG solver_pcg) +FIND_G2O_LIBRARY(G2O_SOLVER_SLAM2D_LINEAR solver_slam2d_linear) +FIND_G2O_LIBRARY(G2O_SOLVER_STRUCTURE_ONLY solver_structure_only) +FIND_G2O_LIBRARY(G2O_SOLVER_EIGEN solver_eigen) + +# Find the predefined types +FIND_G2O_LIBRARY(G2O_TYPES_DATA types_data) +FIND_G2O_LIBRARY(G2O_TYPES_ICP types_icp) +FIND_G2O_LIBRARY(G2O_TYPES_SBA types_sba) +FIND_G2O_LIBRARY(G2O_TYPES_SCLAM2D types_sclam2d) +FIND_G2O_LIBRARY(G2O_TYPES_SIM3 types_sim3) +FIND_G2O_LIBRARY(G2O_TYPES_SLAM2D types_slam2d) +FIND_G2O_LIBRARY(G2O_TYPES_SLAM3D types_slam3d) +FIND_G2O_LIBRARY(G2O_TYPES_SLAM3D_ADDONS types_slam3d_addons) + +# G2O solvers declared found if we found at least one solver +SET(G2O_SOLVERS_FOUND "NO") +IF(G2O_SOLVER_CHOLMOD OR G2O_SOLVER_CSPARSE OR G2O_SOLVER_DENSE OR G2O_SOLVER_PCG OR G2O_SOLVER_SLAM2D_LINEAR OR G2O_SOLVER_STRUCTURE_ONLY OR G2O_SOLVER_EIGEN) + SET(G2O_SOLVERS_FOUND "YES") +ENDIF(G2O_SOLVER_CHOLMOD OR G2O_SOLVER_CSPARSE OR G2O_SOLVER_DENSE OR G2O_SOLVER_PCG OR G2O_SOLVER_SLAM2D_LINEAR OR G2O_SOLVER_STRUCTURE_ONLY OR G2O_SOLVER_EIGEN) + +# G2O itself declared found if we found the core libraries and at least one solver +SET(G2O_FOUND "NO") +IF(G2O_STUFF_LIBRARY AND G2O_CORE_LIBRARY AND G2O_INCLUDE_DIR AND G2O_SOLVERS_FOUND) + SET(G2O_FOUND "YES") +ENDIF(G2O_STUFF_LIBRARY AND G2O_CORE_LIBRARY AND G2O_INCLUDE_DIR AND G2O_SOLVERS_FOUND) diff --git a/docker/build.sh b/docker/build.sh new file mode 100755 index 0000000..17fb7a6 --- /dev/null +++ b/docker/build.sh @@ -0,0 +1,2 @@ +#!/bin/bash +docker build --tag hdl_graph_slam -f noetic/Dockerfile .. diff --git a/docker/howtouse.md b/docker/howtouse.md new file mode 100644 index 0000000..fb5a373 --- /dev/null +++ b/docker/howtouse.md @@ -0,0 +1,40 @@ +# hdl_graph_slam + +Original repository: https://github.com/koide3/hdl_graph_slam + + +## Build +```bash +cd hdl_graph_slam/docker +./build.sh +``` + +## Run + +### On host: +```bash +roscore +``` + +```bash +rosparam set use_sim_time true + +cd hdl_graph_slam/rviz +rviz -d hdl_graph_slam.rviz +``` + +```bash +rosbag play --clock hdl_400.bag +``` +http://www.aisl.cs.tut.ac.jp/databases/hdl_graph_slam/hdl_400.bag.tar.gz + +### On docker image: +```bash +cd hdl_graph_slam/docker +./run.sh + +roslaunch hdl_graph_slam hdl_graph_slam_400.launch +``` + + +![hdl_graph_slam](https://user-images.githubusercontent.com/31344317/98347836-4fed5a00-205b-11eb-931c-158f6cd056bf.gif) diff --git a/docker/kinetic/Dockerfile b/docker/kinetic/Dockerfile new file mode 100644 index 0000000..5331b1c --- /dev/null +++ b/docker/kinetic/Dockerfile @@ -0,0 +1,25 @@ +FROM ros:kinetic + +RUN apt-get update && apt-get install --no-install-recommends -y && apt-get install --no-install-recommends -y wget nano build-essential ros-kinetic-geodesy ros-kinetic-pcl-ros ros-kinetic-nmea-msgs ros-kinetic-rviz ros-kinetic-tf-conversions ros-kinetic-libg2o \ + && apt-get clean \ + && rm -rf /var/lib/apt/lists/* + + +RUN mkdir -p /root/catkin_ws/src +WORKDIR /root/catkin_ws/src +RUN /bin/bash -c '. /opt/ros/kinetic/setup.bash; catkin_init_workspace' +RUN git clone https://github.com/koide3/ndt_omp.git +RUN git clone https://github.com/SMRT-AIST/fast_gicp.git --recursive + +COPY . /root/catkin_ws/src/hdl_graph_slam/ + +WORKDIR /root/catkin_ws +RUN /bin/bash -c '. /opt/ros/kinetic/setup.bash; catkin_make' +RUN sed -i "6i source \"/root/catkin_ws/devel/setup.bash\"" /ros_entrypoint.sh + +RUN apt-get clean && rm -rf /var/lib/apt/lists/* + +WORKDIR / + +ENTRYPOINT ["/ros_entrypoint.sh"] +CMD ["bash"] diff --git a/docker/kinetic_llvm/Dockerfile b/docker/kinetic_llvm/Dockerfile new file mode 100644 index 0000000..e16aa86 --- /dev/null +++ b/docker/kinetic_llvm/Dockerfile @@ -0,0 +1,30 @@ +FROM ros:kinetic + +RUN apt-get update && apt-get install --no-install-recommends -y \ + && apt-get install --no-install-recommends -y \ + wget nano build-essential libomp-dev clang lld-6.0 \ + ros-kinetic-geodesy ros-kinetic-pcl-ros ros-kinetic-nmea-msgs \ + ros-kinetic-rviz ros-kinetic-tf-conversions ros-kinetic-libg2o \ + && apt-get clean \ + && rm -rf /var/lib/apt/lists/* + +RUN update-alternatives --install /usr/bin/ld ld /usr/bin/ld.lld-6.0 50 + +RUN mkdir -p /root/catkin_ws/src +WORKDIR /root/catkin_ws/src +RUN /bin/bash -c '. /opt/ros/kinetic/setup.bash; catkin_init_workspace' +RUN git clone https://github.com/koide3/ndt_omp.git +RUN git clone https://github.com/SMRT-AIST/fast_gicp.git --recursive + +COPY . /root/catkin_ws/src/hdl_graph_slam/ + +WORKDIR /root/catkin_ws +RUN /bin/bash -c '. /opt/ros/kinetic/setup.bash; CC=clang CXX=clang++ catkin_make' +RUN sed -i "6i source \"/root/catkin_ws/devel/setup.bash\"" /ros_entrypoint.sh + +RUN apt-get clean && rm -rf /var/lib/apt/lists/* + +WORKDIR / + +ENTRYPOINT ["/ros_entrypoint.sh"] +CMD ["bash"] diff --git a/docker/melodic/Dockerfile b/docker/melodic/Dockerfile new file mode 100644 index 0000000..0c2237b --- /dev/null +++ b/docker/melodic/Dockerfile @@ -0,0 +1,26 @@ +FROM ros:melodic + +RUN apt-get update && apt-get install -y --no-install-recommends \ + && apt-get install -y --no-install-recommends \ + wget nano build-essential libomp-dev clang lld\ + ros-melodic-geodesy ros-melodic-pcl-ros ros-melodic-nmea-msgs \ + ros-melodic-rviz ros-melodic-tf-conversions ros-melodic-libg2o \ + && apt-get clean \ + && rm -rf /var/lib/apt/lists/* + +RUN mkdir -p /root/catkin_ws/src +WORKDIR /root/catkin_ws/src +RUN /bin/bash -c '. /opt/ros/melodic/setup.bash; catkin_init_workspace' +RUN git clone https://github.com/koide3/ndt_omp.git +RUN git clone https://github.com/SMRT-AIST/fast_gicp.git --recursive + +COPY . /root/catkin_ws/src/hdl_graph_slam/ + +WORKDIR /root/catkin_ws +RUN /bin/bash -c '. /opt/ros/melodic/setup.bash; catkin_make' +RUN sed -i "6i source \"/root/catkin_ws/devel/setup.bash\"" /ros_entrypoint.sh + +WORKDIR / + +ENTRYPOINT ["/ros_entrypoint.sh"] +CMD ["bash"] diff --git a/docker/melodic_llvm/Dockerfile b/docker/melodic_llvm/Dockerfile new file mode 100644 index 0000000..e388b28 --- /dev/null +++ b/docker/melodic_llvm/Dockerfile @@ -0,0 +1,28 @@ +FROM ros:melodic + +RUN apt-get update && apt-get install -y --no-install-recommends \ + && apt-get install -y --no-install-recommends \ + wget nano build-essential libomp-dev clang lld\ + ros-melodic-geodesy ros-melodic-pcl-ros ros-melodic-nmea-msgs \ + ros-melodic-rviz ros-melodic-tf-conversions ros-melodic-libg2o \ + && apt-get clean \ + && rm -rf /var/lib/apt/lists/* + +RUN mkdir -p /root/catkin_ws/src +WORKDIR /root/catkin_ws/src +RUN /bin/bash -c '. /opt/ros/melodic/setup.bash; catkin_init_workspace' +RUN git clone -b melodic https://github.com/koide3/ndt_omp.git +RUN git clone https://github.com/SMRT-AIST/fast_gicp.git --recursive + +COPY . /root/catkin_ws/src/hdl_graph_slam/ + +RUN update-alternatives --install /usr/bin/ld ld /usr/bin/ld.lld 50 + +WORKDIR /root/catkin_ws +RUN /bin/bash -c '. /opt/ros/melodic/setup.bash; CC=clang CXX=clang++ catkin_make' +RUN sed -i "6i source \"/root/catkin_ws/devel/setup.bash\"" /ros_entrypoint.sh + +WORKDIR / + +ENTRYPOINT ["/ros_entrypoint.sh"] +CMD ["bash"] diff --git a/docker/noetic/Dockerfile b/docker/noetic/Dockerfile new file mode 100644 index 0000000..a0759d5 --- /dev/null +++ b/docker/noetic/Dockerfile @@ -0,0 +1,26 @@ +FROM ros:noetic + +RUN apt-get update && apt-get install -y --no-install-recommends \ + && apt-get install -y --no-install-recommends \ + wget nano build-essential libomp-dev clang lld git\ + ros-noetic-geodesy ros-noetic-pcl-ros ros-noetic-nmea-msgs \ + ros-noetic-rviz ros-noetic-tf-conversions ros-noetic-libg2o \ + && apt-get clean \ + && rm -rf /var/lib/apt/lists/* + +RUN mkdir -p /root/catkin_ws/src +WORKDIR /root/catkin_ws/src +RUN /bin/bash -c '. /opt/ros/noetic/setup.bash; catkin_init_workspace' +RUN git clone https://github.com/koide3/ndt_omp.git +RUN git clone https://github.com/SMRT-AIST/fast_gicp.git --recursive + +COPY . /root/catkin_ws/src/hdl_graph_slam/ + +WORKDIR /root/catkin_ws +RUN /bin/bash -c '. /opt/ros/noetic/setup.bash; catkin_make' +RUN sed -i "6i source \"/root/catkin_ws/devel/setup.bash\"" /ros_entrypoint.sh + +WORKDIR / + +ENTRYPOINT ["/ros_entrypoint.sh"] +CMD ["bash"] diff --git a/docker/noetic_llvm/Dockerfile b/docker/noetic_llvm/Dockerfile new file mode 100644 index 0000000..0cfcfbc --- /dev/null +++ b/docker/noetic_llvm/Dockerfile @@ -0,0 +1,28 @@ +FROM ros:noetic + +RUN apt-get update && apt-get install -y --no-install-recommends \ + && apt-get install -y --no-install-recommends \ + wget nano build-essential libomp-dev clang lld git\ + ros-noetic-geodesy ros-noetic-pcl-ros ros-noetic-nmea-msgs \ + ros-noetic-rviz ros-noetic-tf-conversions ros-noetic-libg2o \ + && apt-get clean \ + && rm -rf /var/lib/apt/lists/* + +RUN mkdir -p /root/catkin_ws/src +WORKDIR /root/catkin_ws/src +RUN /bin/bash -c '. /opt/ros/noetic/setup.bash; catkin_init_workspace' +RUN git clone https://github.com/koide3/ndt_omp.git +RUN git clone https://github.com/SMRT-AIST/fast_gicp.git --recursive + +COPY . /root/catkin_ws/src/hdl_graph_slam/ + +RUN update-alternatives --install /usr/bin/ld ld /usr/bin/ld.lld 50 + +WORKDIR /root/catkin_ws +RUN /bin/bash -c '. /opt/ros/noetic/setup.bash; CC=clang CXX=clang++ catkin_make' +RUN sed -i "6i source \"/root/catkin_ws/devel/setup.bash\"" /ros_entrypoint.sh + +WORKDIR / + +ENTRYPOINT ["/ros_entrypoint.sh"] +CMD ["bash"] diff --git a/docker/run.sh b/docker/run.sh new file mode 100755 index 0000000..473cd43 --- /dev/null +++ b/docker/run.sh @@ -0,0 +1,6 @@ +#!/bin/bash +docker run --net=host -it --rm \ + -v $(realpath ..):/root/catkin_ws/src/hdl_graph_slam \ + -w /root/catkin_ws/src/hdl_graph_slam \ + $@ \ + hdl_graph_slam diff --git a/imgs/birds.png b/imgs/birds.png new file mode 100644 index 0000000..8358390 Binary files /dev/null and b/imgs/birds.png differ diff --git a/imgs/ford1.png b/imgs/ford1.png new file mode 100644 index 0000000..b5e4c79 Binary files /dev/null and b/imgs/ford1.png differ diff --git a/imgs/ford2.png b/imgs/ford2.png new file mode 100644 index 0000000..275df01 Binary files /dev/null and b/imgs/ford2.png differ diff --git a/imgs/ford3.png b/imgs/ford3.png new file mode 100644 index 0000000..016b09e Binary files /dev/null and b/imgs/ford3.png differ diff --git a/imgs/hdl_400_graph.png b/imgs/hdl_400_graph.png new file mode 100644 index 0000000..14f5840 Binary files /dev/null and b/imgs/hdl_400_graph.png differ diff --git a/imgs/hdl_400_points.png b/imgs/hdl_400_points.png new file mode 100644 index 0000000..3d08a19 Binary files /dev/null and b/imgs/hdl_400_points.png differ diff --git a/imgs/hdl_graph_slam.png b/imgs/hdl_graph_slam.png new file mode 100644 index 0000000..04c226c Binary files /dev/null and b/imgs/hdl_graph_slam.png differ diff --git a/imgs/nodelets.png b/imgs/nodelets.png new file mode 100644 index 0000000..0d68c2c Binary files /dev/null and b/imgs/nodelets.png differ diff --git a/imgs/nodelets.vsdx b/imgs/nodelets.vsdx new file mode 100644 index 0000000..3311c05 Binary files /dev/null and b/imgs/nodelets.vsdx differ diff --git a/imgs/packages.png b/imgs/packages.png new file mode 100644 index 0000000..b881d70 Binary files /dev/null and b/imgs/packages.png differ diff --git a/imgs/top.png b/imgs/top.png new file mode 100644 index 0000000..df207a8 Binary files /dev/null and b/imgs/top.png differ diff --git a/include/g2o/edge_plane_identity.hpp b/include/g2o/edge_plane_identity.hpp new file mode 100644 index 0000000..b83e4f0 --- /dev/null +++ b/include/g2o/edge_plane_identity.hpp @@ -0,0 +1,103 @@ +// g2o - General Graph Optimization +// Copyright (C) 2011 R. Kuemmerle, G. Grisetti, W. Burgard +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +#ifndef EDGE_PLANE_IDENTITY_HPP +#define EDGE_PLANE_IDENTITY_HPP + +#include +#include +#include + +namespace g2o { + +/** + * @brief A modified version of g2o::EdgePlane. This class takes care of flipped plane normals. + * + */ +class EdgePlaneIdentity : public BaseBinaryEdge<4, Eigen::Vector4d, VertexPlane, VertexPlane> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgePlaneIdentity() : BaseBinaryEdge<4, Eigen::Vector4d, VertexPlane, VertexPlane>() { + _information.setIdentity(); + _error.setZero(); + } + void computeError() { + const VertexPlane* v1 = static_cast(_vertices[0]); + const VertexPlane* v2 = static_cast(_vertices[1]); + + Eigen::Vector4d p1 = v1->estimate().toVector(); + Eigen::Vector4d p2 = v2->estimate().toVector(); + + if(p1.dot(p2) < 0.0) { + p2 = -p2; + } + + _error = (p2 - p1) - _measurement; + } + virtual bool read(std::istream& is) override { + Eigen::Vector4d v; + for(int i = 0; i < 4; ++i) { + is >> v[i]; + } + + setMeasurement(v); + for(int i = 0; i < information().rows(); ++i) { + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) { + information()(j, i) = information()(i, j); + } + } + } + + return true; + } + + virtual bool write(std::ostream& os) const override { + for(int i = 0; i < 4; ++i) { + os << _measurement[i] << " "; + } + + for(int i = 0; i < information().rows(); ++i) { + for(int j = i; j < information().cols(); ++j) { + os << " " << information()(i, j); + }; + } + return os.good(); + } + + virtual void setMeasurement(const Eigen::Vector4d& m) override { + _measurement = m; + } + + virtual int measurementDimension() const override { + return 4; + } +}; + +} // namespace g2o + +#endif // EDGE_PLANE_PARALLEL_HPP diff --git a/include/g2o/edge_plane_parallel.hpp b/include/g2o/edge_plane_parallel.hpp new file mode 100644 index 0000000..edc9ebb --- /dev/null +++ b/include/g2o/edge_plane_parallel.hpp @@ -0,0 +1,159 @@ +// g2o - General Graph Optimization +// Copyright (C) 2011 R. Kuemmerle, G. Grisetti, W. Burgard +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +#ifndef EDGE_PLANE_PARALLEL_HPP +#define EDGE_PLANE_PARALLEL_HPP + +#include +#include +#include + +namespace g2o { + +class EdgePlaneParallel : public BaseBinaryEdge<3, Eigen::Vector3d, VertexPlane, VertexPlane> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgePlaneParallel() : BaseBinaryEdge<3, Eigen::Vector3d, VertexPlane, VertexPlane>() { + _information.setIdentity(); + _error.setZero(); + } + + void computeError() override { + const VertexPlane* v1 = static_cast(_vertices[0]); + const VertexPlane* v2 = static_cast(_vertices[1]); + + Eigen::Vector3d normal1 = v1->estimate().normal(); + Eigen::Vector3d normal2 = v2->estimate().normal(); + + if(normal1.dot(normal2) < 0.0) { + normal2 = -normal2; + } + + _error = (normal2 - normal1) - _measurement; + } + virtual bool read(std::istream& is) override { + Eigen::Vector3d v; + for(int i = 0; i < 3; ++i) { + is >> v[i]; + } + + setMeasurement(v); + for(int i = 0; i < information().rows(); ++i) { + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) { + information()(j, i) = information()(i, j); + } + } + } + + return true; + } + + virtual bool write(std::ostream& os) const override { + for(int i = 0; i < 3; ++i) { + os << _measurement[i] << " "; + } + + for(int i = 0; i < information().rows(); ++i) { + for(int j = i; j < information().cols(); ++j) { + os << " " << information()(i, j); + }; + } + return os.good(); + } + + virtual void setMeasurement(const Eigen::Vector3d& m) override { + _measurement = m; + } + + virtual int measurementDimension() const override { + return 3; + } +}; + +class EdgePlanePerpendicular : public BaseBinaryEdge<1, Eigen::Vector3d, VertexPlane, VertexPlane> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgePlanePerpendicular() : BaseBinaryEdge<1, Eigen::Vector3d, VertexPlane, VertexPlane>() { + _information.setIdentity(); + _error.setZero(); + } + + void computeError() override { + const VertexPlane* v1 = static_cast(_vertices[0]); + const VertexPlane* v2 = static_cast(_vertices[1]); + + Eigen::Vector3d normal1 = v1->estimate().normal().normalized(); + Eigen::Vector3d normal2 = v2->estimate().normal().normalized(); + + _error[0] = normal1.dot(normal2); + } + + virtual bool read(std::istream& is) override { + Eigen::Vector3d v; + for(int i = 0; i < 3; ++i) { + is >> v[i]; + } + + setMeasurement(v); + for(int i = 0; i < information().rows(); ++i) { + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) { + information()(j, i) = information()(i, j); + } + } + } + + return true; + } + + virtual bool write(std::ostream& os) const override { + for(int i = 0; i < 3; ++i) { + os << _measurement[i] << " "; + } + + for(int i = 0; i < information().rows(); ++i) { + for(int j = i; j < information().cols(); ++j) { + os << " " << information()(i, j); + }; + } + return os.good(); + } + + virtual void setMeasurement(const Eigen::Vector3d& m) override { + _measurement = m; + } + + virtual int measurementDimension() const override { + return 3; + } +}; + +} // namespace g2o + +#endif // EDGE_PLANE_PARALLEL_HPP diff --git a/include/g2o/edge_plane_prior.hpp b/include/g2o/edge_plane_prior.hpp new file mode 100644 index 0000000..cd914af --- /dev/null +++ b/include/g2o/edge_plane_prior.hpp @@ -0,0 +1,107 @@ +// g2o - General Graph Optimization +// Copyright (C) 2011 R. Kuemmerle, G. Grisetti, W. Burgard +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +#ifndef EDGE_PLANE_PRIOR_HPP +#define EDGE_PLANE_PRIOR_HPP + +#include +#include +#include + +namespace g2o { +class EdgePlanePriorNormal : public g2o::BaseUnaryEdge<3, Eigen::Vector3d, g2o::VertexPlane> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgePlanePriorNormal() : g2o::BaseUnaryEdge<3, Eigen::Vector3d, g2o::VertexPlane>() {} + + void computeError() override { + const g2o::VertexPlane* v1 = static_cast(_vertices[0]); + Eigen::Vector3d normal = v1->estimate().normal(); + + if(normal.dot(_measurement) < 0.0) { + normal = -normal; + } + + _error = normal - _measurement; + } + + void setMeasurement(const Eigen::Vector3d& m) override { + _measurement = m; + } + + virtual bool read(std::istream& is) override { + Eigen::Vector3d v; + is >> v(0) >> v(1) >> v(2); + setMeasurement(v); + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) information()(j, i) = information()(i, j); + } + return true; + } + virtual bool write(std::ostream& os) const override { + Eigen::Vector3d v = _measurement; + os << v(0) << " " << v(1) << " " << v(2) << " "; + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) os << " " << information()(i, j); + return os.good(); + } +}; + +class EdgePlanePriorDistance : public g2o::BaseUnaryEdge<1, double, g2o::VertexPlane> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgePlanePriorDistance() : g2o::BaseUnaryEdge<1, double, g2o::VertexPlane>() {} + + void computeError() override { + const g2o::VertexPlane* v1 = static_cast(_vertices[0]); + _error[0] = _measurement - v1->estimate().distance(); + } + + void setMeasurement(const double& m) override { + _measurement = m; + } + + virtual bool read(std::istream& is) override { + is >> _measurement; + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) information()(j, i) = information()(i, j); + } + return true; + } + virtual bool write(std::ostream& os) const override { + os << _measurement; + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) os << " " << information()(i, j); + return os.good(); + } +}; +} // namespace g2o + +#endif // EDGE_SE3_PRIORXY_HPP diff --git a/include/g2o/edge_se3_plane.hpp b/include/g2o/edge_se3_plane.hpp new file mode 100644 index 0000000..bf3dedb --- /dev/null +++ b/include/g2o/edge_se3_plane.hpp @@ -0,0 +1,74 @@ +// g2o - General Graph Optimization +// Copyright (C) 2011 R. Kuemmerle, G. Grisetti, W. Burgard +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +#ifndef KKL_G2O_EDGE_SE3_PLANE_HPP +#define KKL_G2O_EDGE_SE3_PLANE_HPP + +#include +#include +#include + +namespace g2o { +class EdgeSE3Plane : public g2o::BaseBinaryEdge<3, g2o::Plane3D, g2o::VertexSE3, g2o::VertexPlane> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgeSE3Plane() : BaseBinaryEdge<3, g2o::Plane3D, g2o::VertexSE3, g2o::VertexPlane>() {} + + void computeError() override { + const g2o::VertexSE3* v1 = static_cast(_vertices[0]); + const g2o::VertexPlane* v2 = static_cast(_vertices[1]); + + Eigen::Isometry3d w2n = v1->estimate().inverse(); + Plane3D local_plane = w2n * v2->estimate(); + _error = local_plane.ominus(_measurement); + } + + void setMeasurement(const g2o::Plane3D& m) override { + _measurement = m; + } + + virtual bool read(std::istream& is) override { + Eigen::Vector4d v; + is >> v(0) >> v(1) >> v(2) >> v(3); + setMeasurement(Plane3D(v)); + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) information()(j, i) = information()(i, j); + } + return true; + } + virtual bool write(std::ostream& os) const override { + Eigen::Vector4d v = _measurement.toVector(); + os << v(0) << " " << v(1) << " " << v(2) << " " << v(3) << " "; + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) os << " " << information()(i, j); + return os.good(); + } +}; +} // namespace g2o + +#endif diff --git a/include/g2o/edge_se3_priorquat.hpp b/include/g2o/edge_se3_priorquat.hpp new file mode 100644 index 0000000..f05f34e --- /dev/null +++ b/include/g2o/edge_se3_priorquat.hpp @@ -0,0 +1,78 @@ +// g2o - General Graph Optimization +// Copyright (C) 2011 R. Kuemmerle, G. Grisetti, W. Burgard +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +#ifndef KKL_G2O_EDGE_SE3_PRIORQUAT_HPP +#define KKL_G2O_EDGE_SE3_PRIORQUAT_HPP + +#include +#include + +namespace g2o { +class EdgeSE3PriorQuat : public g2o::BaseUnaryEdge<3, Eigen::Quaterniond, g2o::VertexSE3> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgeSE3PriorQuat() : g2o::BaseUnaryEdge<3, Eigen::Quaterniond, g2o::VertexSE3>() {} + + void computeError() override { + const g2o::VertexSE3* v1 = static_cast(_vertices[0]); + + Eigen::Quaterniond estimate = Eigen::Quaterniond(v1->estimate().linear()); + + if(_measurement.coeffs().dot(estimate.coeffs()) < 0.0) { + estimate.coeffs() = -estimate.coeffs(); + } + _error = estimate.vec() - _measurement.vec(); + } + + void setMeasurement(const Eigen::Quaterniond& m) override { + _measurement = m; + if(m.w() < 0.0) { + _measurement.coeffs() = -m.coeffs(); + } + } + + virtual bool read(std::istream& is) override { + Eigen::Quaterniond q; + is >> q.w() >> q.x() >> q.y() >> q.z(); + setMeasurement(q); + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) information()(j, i) = information()(i, j); + } + return true; + } + virtual bool write(std::ostream& os) const override { + Eigen::Quaterniond q = _measurement; + os << q.w() << " " << q.x() << " " << q.y() << " " << q.z(); + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) os << " " << information()(i, j); + return os.good(); + } +}; +} // namespace g2o + +#endif diff --git a/include/g2o/edge_se3_priorvec.hpp b/include/g2o/edge_se3_priorvec.hpp new file mode 100644 index 0000000..4d2115f --- /dev/null +++ b/include/g2o/edge_se3_priorvec.hpp @@ -0,0 +1,76 @@ +// g2o - General Graph Optimization +// Copyright (C) 2011 R. Kuemmerle, G. Grisetti, W. Burgard +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +#ifndef KKL_G2O_EDGE_SE3_PRIORVEC_HPP +#define KKL_G2O_EDGE_SE3_PRIORVEC_HPP + +#include +#include + +namespace g2o { +class EdgeSE3PriorVec : public g2o::BaseUnaryEdge<3, Eigen::Matrix, g2o::VertexSE3> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgeSE3PriorVec() : g2o::BaseUnaryEdge<3, Eigen::Matrix, g2o::VertexSE3>() {} + + void computeError() override { + const g2o::VertexSE3* v1 = static_cast(_vertices[0]); + + Eigen::Vector3d direction = _measurement.head<3>(); + Eigen::Vector3d measurement = _measurement.tail<3>(); + + Eigen::Vector3d estimate = (v1->estimate().linear().inverse() * direction); + + _error = estimate - measurement; + } + + void setMeasurement(const Eigen::Matrix& m) override { + _measurement.head<3>() = m.head<3>().normalized(); + _measurement.tail<3>() = m.tail<3>().normalized(); + } + + virtual bool read(std::istream& is) override { + Eigen::Matrix v; + is >> v[0] >> v[1] >> v[2] >> v[3] >> v[4] >> v[5]; + setMeasurement(v); + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) information()(j, i) = information()(i, j); + } + return true; + } + virtual bool write(std::ostream& os) const override { + Eigen::Matrix v = _measurement; + os << v[0] << " " << v[1] << " " << v[2] << " " << v[3] << " " << v[4] << " " << v[5]; + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) os << " " << information()(i, j); + return os.good(); + } +}; +} // namespace g2o + +#endif diff --git a/include/g2o/edge_se3_priorxy.hpp b/include/g2o/edge_se3_priorxy.hpp new file mode 100644 index 0000000..d8fb35c --- /dev/null +++ b/include/g2o/edge_se3_priorxy.hpp @@ -0,0 +1,71 @@ +// g2o - General Graph Optimization +// Copyright (C) 2011 R. Kuemmerle, G. Grisetti, W. Burgard +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +#ifndef EDGE_SE3_PRIORXY_HPP +#define EDGE_SE3_PRIORXY_HPP + +#include +#include + +namespace g2o { +class EdgeSE3PriorXY : public g2o::BaseUnaryEdge<2, Eigen::Vector2d, g2o::VertexSE3> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgeSE3PriorXY() : g2o::BaseUnaryEdge<2, Eigen::Vector2d, g2o::VertexSE3>() {} + + void computeError() override { + const g2o::VertexSE3* v1 = static_cast(_vertices[0]); + + Eigen::Vector2d estimate = v1->estimate().translation().head<2>(); + _error = estimate - _measurement; + } + + void setMeasurement(const Eigen::Vector2d& m) override { + _measurement = m; + } + + virtual bool read(std::istream& is) override { + Eigen::Vector2d v; + is >> v(0) >> v(1); + setMeasurement(v); + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) information()(j, i) = information()(i, j); + } + return true; + } + virtual bool write(std::ostream& os) const override { + Eigen::Vector2d v = _measurement; + os << v(0) << " " << v(1) << " "; + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) os << " " << information()(i, j); + return os.good(); + } +}; +} // namespace g2o + +#endif // EDGE_SE3_PRIORXY_HPP diff --git a/include/g2o/edge_se3_priorxyz.hpp b/include/g2o/edge_se3_priorxyz.hpp new file mode 100644 index 0000000..575f924 --- /dev/null +++ b/include/g2o/edge_se3_priorxyz.hpp @@ -0,0 +1,71 @@ +// g2o - General Graph Optimization +// Copyright (C) 2011 R. Kuemmerle, G. Grisetti, W. Burgard +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +// IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED +// TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +#ifndef KKL_G2O_EDGE_SE3_PRIORXYZ_HPP +#define KKL_G2O_EDGE_SE3_PRIORXYZ_HPP + +#include +#include + +namespace g2o { +class EdgeSE3PriorXYZ : public g2o::BaseUnaryEdge<3, Eigen::Vector3d, g2o::VertexSE3> { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + EdgeSE3PriorXYZ() : g2o::BaseUnaryEdge<3, Eigen::Vector3d, g2o::VertexSE3>() {} + + void computeError() override { + const g2o::VertexSE3* v1 = static_cast(_vertices[0]); + + Eigen::Vector3d estimate = v1->estimate().translation(); + _error = estimate - _measurement; + } + + void setMeasurement(const Eigen::Vector3d& m) override { + _measurement = m; + } + + virtual bool read(std::istream& is) override { + Eigen::Vector3d v; + is >> v(0) >> v(1) >> v(2); + setMeasurement(v); + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) { + is >> information()(i, j); + if(i != j) information()(j, i) = information()(i, j); + } + return true; + } + virtual bool write(std::ostream& os) const override { + Eigen::Vector3d v = _measurement; + os << v(0) << " " << v(1) << " " << v(2) << " "; + for(int i = 0; i < information().rows(); ++i) + for(int j = i; j < information().cols(); ++j) os << " " << information()(i, j); + return os.good(); + } +}; +} // namespace g2o + +#endif diff --git a/include/g2o/robust_kernel_io.hpp b/include/g2o/robust_kernel_io.hpp new file mode 100644 index 0000000..ed5a53b --- /dev/null +++ b/include/g2o/robust_kernel_io.hpp @@ -0,0 +1,19 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef ROBUST_KERNEL_IO_HPP +#define ROBUST_KERNEL_IO_HPP + +#include +#include + +namespace g2o { + +std::string kernel_type(g2o::RobustKernel* kernel); + +bool save_robust_kernels(const std::string& filename, g2o::SparseOptimizer* graph); + +bool load_robust_kernels(const std::string& filename, g2o::SparseOptimizer* graph); + +} // namespace g2o + +#endif // ROBUST_KERNEL_IO_HPP diff --git a/include/hdl_graph_slam/featureAssociation.h b/include/hdl_graph_slam/featureAssociation.h new file mode 100644 index 0000000..6774bd2 --- /dev/null +++ b/include/hdl_graph_slam/featureAssociation.h @@ -0,0 +1,236 @@ +// +// Created by robot410 on 2020/11/3. +// + +#ifndef LOAM_ODOMETRY_FEATUREASSOCIATION_H +#define LOAM_ODOMETRY_FEATUREASSOCIATION_H + +#include "hdl_graph_slam/utility.h" +typedef Eigen::Matrix Vector21f; +namespace hdl_graph_slam { + class FeatureAssociation { + + public: + + pcl::PointCloud::Ptr segmentedCloud; + pcl::PointCloud::Ptr outlierCloud; + pcl::PointCloud LaserCloud_; + pcl::PointCloud LaserCloud_last_; + + pcl::PointCloud LinearCloud_; + pcl::PointCloud LinearCloud_last_; + + pcl::PointCloud::Ptr Source_filtered_; + pcl::PointCloud::Ptr Target_filtered_; + + pcl::PointCloud Less_LinearCloud_; + pcl::PointCloud Less_LinearCloud_last_; + + pcl::PointCloud LinearCloud_Ground_; + pcl::PointCloud LinearCloud_Less_Ground_; + + pcl::PointCloud LinearCloud_Ground_Last; + pcl::PointCloud LinearCloud_Less_Ground_Last; + + pcl::PointCloud::Ptr cornerPointsSharp; + pcl::PointCloud::Ptr cornerPointsLessSharp; + pcl::PointCloud::Ptr surfPointsFlat; + pcl::PointCloud::Ptr surfPointsLessFlat; + + pcl::PointCloud::Ptr surfPointsLessFlatScan; + pcl::PointCloud::Ptr surfPointsLessFlatScanDS; + + pcl::VoxelGrid downSizeFilter; + + Eigen::Matrix4d RL_transfomation_; + + Eigen::MatrixXd Pose_COV_; + + double timeScanCur; + double timeNewSegmentedCloud; + double timeNewSegmentedCloudInfo; + double timeNewOutlierCloud; + + bool newSegmentedCloud; + bool newSegmentedCloudInfo; + bool newOutlierCloud; + double Time_; + double Thre_; + bool use_sampling_thre_; + + cloud_info segInfo; + + int systemInitCount; + bool systemInited; + + std::vector cloudSmoothness; + float cloudCurvature[N_SCAN * Horizon_SCAN]; + int cloudNeighborPicked[N_SCAN * Horizon_SCAN]; + int cloudLabel[N_SCAN * Horizon_SCAN]; + + int imuPointerFront; + int imuPointerLast; + int imuPointerLastIteration; + + float imuRollStart, imuPitchStart, imuYawStart; + float cosImuRollStart, cosImuPitchStart, cosImuYawStart, sinImuRollStart, sinImuPitchStart, sinImuYawStart; + float imuRollCur, imuPitchCur, imuYawCur; + + float imuVeloXStart, imuVeloYStart, imuVeloZStart; + float imuShiftXStart, imuShiftYStart, imuShiftZStart; + + float imuVeloXCur, imuVeloYCur, imuVeloZCur; + float imuShiftXCur, imuShiftYCur, imuShiftZCur; + + float imuShiftFromStartXCur, imuShiftFromStartYCur, imuShiftFromStartZCur; + float imuVeloFromStartXCur, imuVeloFromStartYCur, imuVeloFromStartZCur; + + float imuAngularRotationXCur, imuAngularRotationYCur, imuAngularRotationZCur; + float imuAngularRotationXLast, imuAngularRotationYLast, imuAngularRotationZLast; + float imuAngularFromStartX, imuAngularFromStartY, imuAngularFromStartZ; + + double imuTime[imuQueLength]; + float imuRoll[imuQueLength]; + float imuPitch[imuQueLength]; + float imuYaw[imuQueLength]; + + float imuAccX[imuQueLength]; + float imuAccY[imuQueLength]; + float imuAccZ[imuQueLength]; + + float imuVeloX[imuQueLength]; + float imuVeloY[imuQueLength]; + float imuVeloZ[imuQueLength]; + + float imuShiftX[imuQueLength]; + float imuShiftY[imuQueLength]; + float imuShiftZ[imuQueLength]; + + float imuAngularVeloX[imuQueLength]; + float imuAngularVeloY[imuQueLength]; + float imuAngularVeloZ[imuQueLength]; + + float imuAngularRotationX[imuQueLength]; + float imuAngularRotationY[imuQueLength]; + float imuAngularRotationZ[imuQueLength]; + + int skipFrameNum; + bool systemInitedLM; + + int laserCloudCornerLastNum; + int laserCloudSurfLastNum; + + int pointSelCornerInd[N_SCAN * Horizon_SCAN]; + float pointSearchCornerInd1[N_SCAN * Horizon_SCAN]; + float pointSearchCornerInd2[N_SCAN * Horizon_SCAN]; + + int pointSelSurfInd[N_SCAN * Horizon_SCAN]; + float pointSearchSurfInd1[N_SCAN * Horizon_SCAN]; + float pointSearchSurfInd2[N_SCAN * Horizon_SCAN]; + float pointSearchSurfInd3[N_SCAN * Horizon_SCAN]; + + float transformCur[6]; + float transformCur_GT[6]; + float transformSum[6]; + float transformIncre[0]; + + float mean_matches_errors_; + float matched_f_no_; + float total_f_no_; + + float Loam_mean_matches_errors_; + float Loam_matched_f_no_; + float Loam_total_f_no_; + + + float imuRollLast, imuPitchLast, imuYawLast; + float imuShiftFromStartX, imuShiftFromStartY, imuShiftFromStartZ; + float imuVeloFromStartX, imuVeloFromStartY, imuVeloFromStartZ; + + pcl::PointCloud::Ptr laserCloudCornerLast; + pcl::PointCloud::Ptr laserCloudSurfLast; + pcl::PointCloud::Ptr laserCloudOri; + pcl::PointCloud::Ptr laserCloudOri_sub; + pcl::PointCloud::Ptr coeffSel; + pcl::PointCloud::Ptr target_cloud_; + + pcl::KdTreeFLANN::Ptr kdtreeLineCornerLast; + pcl::KdTreeFLANN::Ptr kdtreeGroundCornerLast; + pcl::KdTreeFLANN::Ptr kdtreeCornerLast; + pcl::KdTreeFLANN::Ptr kdtreeSurfLast; + + std::vector pointSearchInd; + std::vector pointSearchSqDis; + + vector init_error_; + float a1, a2, a3, a4, a5, a6, a7; + int n1, n2, n3, n4, n5, n6, n7; + vector v1, v2, v3, v4, v5, v6, v7; + vector sigma_set_; + double error_sum; + float avg_smooth_; + float sigma_smooth_; + float sigma_smooth_vs_; + int avg_smooth_no_; + int segmentedCloudNum_; + float avg_error_last_; + Vector21f smooth_dist, smooth_dist_no_; + vector truth_smooth_set_; + vector false_smooth_set_; + float avg_smooth_vs_; + int avg_smooth_no_vs_; + + + float max_prob_; + + double total_cost_ = 0; + Eigen::Matrix4d GT_; + + PointType pointOri, pointSel, tripod1, tripod2, tripod3, pointProj, coeff; + + bool isDegenerate; + cv::Mat matP; + + int iterCount1_, iterCount2_; + int frameCount; + + FeatureAssociation(); + + void TransformToStart(PointType const *const pi, PointType *const po); + + double rad2deg(double radians); + + //double deg2rad(double degrees); + + void findCorrespondingCornerFeatures_X(int iterCount); + + void initializationValue(); + + void updateTransformation(); + + void findCorrespondingSurfFeatures(int iterCount); + + bool calculateTransformationSurf(int iterCount); + + bool calculateTransformationCorner(int iterCount); + + bool calculateTransformation(int iterCount); + + void runFeatureAssociation(Eigen::Matrix4f Init_state); + + void Matrix4dToRpyxyz(const Eigen::Matrix4f &Trans,float nums[]); + + void toEulerAngle( Eigen::Quaternionf& q, float& roll, float& pitch, float& yaw); + + void setSourceCloud(const pcl::PointCloud::ConstPtr cloud); + + void setTargetCloud(const pcl::PointCloud::ConstPtr cloud); + void setInitTargetCloud(const pcl::PointCloud::ConstPtr cloud); + + + Eigen::Matrix4f getFinalTransformation(); + + void RpyxyzToMatrix4f(Eigen::Matrix4f &Trans, float nums[]); + }; +} +#endif //LOAM_ODOMETRY_FEATUREASSOCIATION_H diff --git a/include/hdl_graph_slam/graph_slam.hpp b/include/hdl_graph_slam/graph_slam.hpp new file mode 100644 index 0000000..51efafb --- /dev/null +++ b/include/hdl_graph_slam/graph_slam.hpp @@ -0,0 +1,147 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef GRAPH_SLAM_HPP +#define GRAPH_SLAM_HPP + +#include +#include + +#include + +namespace g2o { +class VertexSE3; +class VertexPlane; +class VertexPointXYZ; +class EdgeSE3; +class EdgeSE3Plane; +class EdgeSE3PointXYZ; +class EdgeSE3PriorXY; +class EdgeSE3PriorXYZ; +class EdgeSE3PriorVec; +class EdgeSE3PriorQuat; +class EdgePlane; +class EdgePlaneIdentity; +class EdgePlaneParallel; +class EdgePlanePerpendicular; +class EdgePlanePriorNormal; +class EdgePlanePriorDistance; +class RobustKernelFactory; +} // namespace g2o + +namespace hdl_graph_slam { + +class GraphSLAM { +public: + GraphSLAM(const std::string& solver_type = "lm_var"); + virtual ~GraphSLAM(); + + int num_vertices() const; + int num_edges() const; + + void set_solver(const std::string& solver_type); + + /** + * @brief add a SE3 node to the graph + * @param pose + * @return registered node + */ + g2o::VertexSE3* add_se3_node(const Eigen::Isometry3d& pose); + + /** + * @brief add a plane node to the graph + * @param plane_coeffs + * @return registered node + */ + g2o::VertexPlane* add_plane_node(const Eigen::Vector4d& plane_coeffs); + + /** + * @brief add a point_xyz node to the graph + * @param xyz + * @return registered node + */ + g2o::VertexPointXYZ* add_point_xyz_node(const Eigen::Vector3d& xyz); + + /** + * @brief add an edge between SE3 nodes + * @param v1 node1 + * @param v2 node2 + * @param relative_pose relative pose between node1 and node2 + * @param information_matrix information matrix (it must be 6x6) + * @return registered edge + */ + g2o::EdgeSE3* add_se3_edge(g2o::VertexSE3* v1, g2o::VertexSE3* v2, const Eigen::Isometry3d& relative_pose, const Eigen::MatrixXd& information_matrix); + + /** + * @brief add an edge between an SE3 node and a plane node + * @param v_se3 SE3 node + * @param v_plane plane node + * @param plane_coeffs plane coefficients w.r.t. v_se3 + * @param information_matrix information matrix (it must be 3x3) + * @return registered edge + */ + g2o::EdgeSE3Plane* add_se3_plane_edge(g2o::VertexSE3* v_se3, g2o::VertexPlane* v_plane, const Eigen::Vector4d& plane_coeffs, const Eigen::MatrixXd& information_matrix); + + /** + * @brief add an edge between an SE3 node and a point_xyz node + * @param v_se3 SE3 node + * @param v_xyz point_xyz node + * @param xyz xyz coordinate + * @param information information_matrix (it must be 3x3) + * @return registered edge + */ + g2o::EdgeSE3PointXYZ* add_se3_point_xyz_edge(g2o::VertexSE3* v_se3, g2o::VertexPointXYZ* v_xyz, const Eigen::Vector3d& xyz, const Eigen::MatrixXd& information_matrix); + + /** + * @brief add a prior edge to an SE3 node + * @param v_se3 + * @param xy + * @param information_matrix + * @return + */ + g2o::EdgePlanePriorNormal* add_plane_normal_prior_edge(g2o::VertexPlane* v, const Eigen::Vector3d& normal, const Eigen::MatrixXd& information_matrix); + + g2o::EdgePlanePriorDistance* add_plane_distance_prior_edge(g2o::VertexPlane* v, double distance, const Eigen::MatrixXd& information_matrix); + + g2o::EdgeSE3PriorXY* add_se3_prior_xy_edge(g2o::VertexSE3* v_se3, const Eigen::Vector2d& xy, const Eigen::MatrixXd& information_matrix); + + g2o::EdgeSE3PriorXYZ* add_se3_prior_xyz_edge(g2o::VertexSE3* v_se3, const Eigen::Vector3d& xyz, const Eigen::MatrixXd& information_matrix); + + g2o::EdgeSE3PriorQuat* add_se3_prior_quat_edge(g2o::VertexSE3* v_se3, const Eigen::Quaterniond& quat, const Eigen::MatrixXd& information_matrix); + + g2o::EdgeSE3PriorVec* add_se3_prior_vec_edge(g2o::VertexSE3* v_se3, const Eigen::Vector3d& direction, const Eigen::Vector3d& measurement, const Eigen::MatrixXd& information_matrix); + + g2o::EdgePlane* add_plane_edge(g2o::VertexPlane* v_plane1, g2o::VertexPlane* v_plane2, const Eigen::Vector4d& measurement, const Eigen::Matrix4d& information); + + g2o::EdgePlaneIdentity* add_plane_identity_edge(g2o::VertexPlane* v_plane1, g2o::VertexPlane* v_plane2, const Eigen::Vector4d& measurement, const Eigen::Matrix4d& information); + + g2o::EdgePlaneParallel* add_plane_parallel_edge(g2o::VertexPlane* v_plane1, g2o::VertexPlane* v_plane2, const Eigen::Vector3d& measurement, const Eigen::Matrix3d& information); + + g2o::EdgePlanePerpendicular* add_plane_perpendicular_edge(g2o::VertexPlane* v_plane1, g2o::VertexPlane* v_plane2, const Eigen::Vector3d& measurement, const Eigen::MatrixXd& information); + + void add_robust_kernel(g2o::HyperGraph::Edge* edge, const std::string& kernel_type, double kernel_size); + + /** + * @brief perform graph optimization + */ + int optimize(int num_iterations); + + /** + * @brief save the pose graph to a file + * @param filename output filename + */ + void save(const std::string& filename); + + /** + * @brief load the pose graph from file + * @param filename output filename + */ + bool load(const std::string& filename); + +public: + g2o::RobustKernelFactory* robust_kernel_factory; + std::unique_ptr graph; // g2o graph +}; + +} // namespace hdl_graph_slam + +#endif // GRAPH_SLAM_HPP diff --git a/include/hdl_graph_slam/imageProjection.h b/include/hdl_graph_slam/imageProjection.h new file mode 100644 index 0000000..68e77bf --- /dev/null +++ b/include/hdl_graph_slam/imageProjection.h @@ -0,0 +1,113 @@ +// +// Created by robot410 on 2020/11/3. +// + +#ifndef LOAM_ODOMETRY_IMAGEPROJECTION_H +#define LOAM_ODOMETRY_IMAGEPROJECTION_H + +#include "hdl_graph_slam/utility.h" +#include +namespace hdl_graph_slam { + typedef Eigen::Matrix Vec16f; + typedef Eigen::Matrix Vector21f; + + class ImageProjection { + public: + + pcl::PointCloud::Ptr laserCloudIn; + pcl::PointCloud::Ptr laserCloudIn_last_; + + pcl::PointCloud::Ptr fullCloud; // projected velodyne raw cloud, but saved in the form of 1-D matrix + pcl::PointCloud::Ptr fullInfoCloud; // same as fullCloud, but with intensity - range + pcl::PointCloud::Ptr LinearPoints_; + pcl::PointCloud::Ptr Less_LinearPoints_; + pcl::PointCloud::Ptr Less_LinearPoints_D_; + pcl::PointCloud::Ptr LinearPoints_last_; + pcl::PointCloud::Ptr fullLinearPoints_; + + pcl::PointCloud::Ptr fullValidPoints_z_; + pcl::PointCloud::Ptr fullValidPoints_smoothness_; + + pcl::PointCloud::Ptr LinearPoints_Ground; + pcl::PointCloud::Ptr LinearPoints_Less_Ground; + + + std::vector pointSearchInd; + std::vector pointSearchSqDis; + + float avg_smoothness_; + float avg_smooth_predict_; + float sigma_smooth_predict_; + vector predicted_sigma_set_; + + float avg_smooth_predict_vs_; + float sigma_smooth_predict_vs_; + + float traget_distribution(float smooth, double s); + + vector neighbor_set_; + + pcl::PointCloud::Ptr groundCloud; + pcl::PointCloud::Ptr segmentedCloud; + pcl::PointCloud::Ptr segmentedCloudPure; + pcl::PointCloud::Ptr outlierCloud; + + pcl::KdTreeFLANN::Ptr kdtree_linear_; + + bool use_prob_smooth_; + + int feature_no_; + + PointType nanPoint; // fill in fullCloud at each iteration + + cv::Mat rangeMat; // range matrix for range image + cv::Mat labelMat; // label matrix for segmentaiton marking + cv::Mat groundMat; // ground matrix for ground cloud marking + int labelCount; + + float startOrientation; + float endOrientation; + + Vector21f last_distr_; + + cloud_info segMsg; // info of segmented cloud + + std::vector > neighborIterator; // neighbor iterator for segmentaiton process + + uint16_t *allPushedIndX; // array for tracking points of a segmented object + uint16_t *allPushedIndY; + + uint16_t *queueIndX; // array for breadth-first search process of segmentation + uint16_t *queueIndY; + + Vec16f smooth_thre_; + + float plane_no_; + + double Time_; + double Thre_; + + pcl::KdTreeFLANN::Ptr kdtree_; + + ImageProjection(); + + Eigen::Matrix4d Init_Trans_; + + void allocateMemory(); + + void resetParameters(); + + void cloudHandler(pcl::PointCloud::Ptr input_cloud); + + void findStartEndAngle(); + + void projectPointCloud(); + + void groundRemoval(); + + void extractLinearPoints(); + + void showTwoPC(pcl::PointCloud::Ptr curr_cloud_gicp, pcl::PointCloud::Ptr ref_cloud); + }; +} +#endif //LOAM_ODOMETRY_IMAGEPROJECTION_H diff --git a/include/hdl_graph_slam/information_matrix_calculator.hpp b/include/hdl_graph_slam/information_matrix_calculator.hpp new file mode 100644 index 0000000..8557fbf --- /dev/null +++ b/include/hdl_graph_slam/information_matrix_calculator.hpp @@ -0,0 +1,59 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef INFORMATION_MATRIX_CALCULATOR_HPP +#define INFORMATION_MATRIX_CALCULATOR_HPP + +#include +#include +#include + +namespace hdl_graph_slam { + +class InformationMatrixCalculator { +public: + using PointT = pcl::PointXYZI; + + InformationMatrixCalculator() {} + InformationMatrixCalculator(ros::NodeHandle& nh); + ~InformationMatrixCalculator(); + + template + void load(ParamServer& params) { + use_const_inf_matrix = params.template param("use_const_inf_matrix", false); + const_stddev_x = params.template param("const_stddev_x", 0.5); + const_stddev_q = params.template param("const_stddev_q", 0.1); + + var_gain_a = params.template param("var_gain_a", 20.0); + min_stddev_x = params.template param("min_stddev_x", 0.1); + max_stddev_x = params.template param("max_stddev_x", 5.0); + min_stddev_q = params.template param("min_stddev_q", 0.05); + max_stddev_q = params.template param("max_stddev_q", 0.2); + fitness_score_thresh = params.template param("fitness_score_thresh", 2.5); + } + + static double calc_fitness_score(const pcl::PointCloud::ConstPtr& cloud1, const pcl::PointCloud::ConstPtr& cloud2, const Eigen::Isometry3d& relpose, double max_range = std::numeric_limits::max()); + + Eigen::MatrixXd calc_information_matrix(const pcl::PointCloud::ConstPtr& cloud1, const pcl::PointCloud::ConstPtr& cloud2, const Eigen::Isometry3d& relpose) const; + +private: + double weight(double a, double max_x, double min_y, double max_y, double x) const { + double y = (1.0 - std::exp(-a * x)) / (1.0 - std::exp(-a * max_x)); + return min_y + (max_y - min_y) * y; + } + +private: + bool use_const_inf_matrix; + double const_stddev_x; + double const_stddev_q; + + double var_gain_a; + double min_stddev_x; + double max_stddev_x; + double min_stddev_q; + double max_stddev_q; + double fitness_score_thresh; +}; + +} // namespace hdl_graph_slam + +#endif // INFORMATION_MATRIX_CALCULATOR_HPP diff --git a/include/hdl_graph_slam/keyframe.hpp b/include/hdl_graph_slam/keyframe.hpp new file mode 100644 index 0000000..8714d43 --- /dev/null +++ b/include/hdl_graph_slam/keyframe.hpp @@ -0,0 +1,74 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef KEYFRAME_HPP +#define KEYFRAME_HPP + +#include +#include +#include +#include + +namespace g2o { +class VertexSE3; +class HyperGraph; +class SparseOptimizer; +} // namespace g2o + +namespace hdl_graph_slam { + +/** + * @brief KeyFrame (pose node) + */ +struct KeyFrame { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + using PointT = pcl::PointXYZI; + using Ptr = std::shared_ptr; + + KeyFrame(const ros::Time& stamp, const Eigen::Isometry3d& odom, double accum_distance, const pcl::PointCloud::ConstPtr& cloud); + KeyFrame(const std::string& directory, g2o::HyperGraph* graph); + virtual ~KeyFrame(); + + void save(const std::string& directory); + bool load(const std::string& directory, g2o::HyperGraph* graph); + + long id() const; + Eigen::Isometry3d estimate() const; + +public: + ros::Time stamp; // timestamp + Eigen::Isometry3d odom; // odometry (estimated by scan_matching_odometry) + double accum_distance; // accumulated distance from the first node (by scan_matching_odometry) + pcl::PointCloud::ConstPtr cloud; // point cloud + boost::optional floor_coeffs; // detected floor's coefficients + boost::optional utm_coord; // UTM coord obtained by GPS + + boost::optional acceleration; // + boost::optional orientation; // + + g2o::VertexSE3* node; // node instance +}; + +/** + * @brief KeyFramesnapshot for map cloud generation + */ +struct KeyFrameSnapshot { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + + using PointT = KeyFrame::PointT; + using Ptr = std::shared_ptr; + + KeyFrameSnapshot(const KeyFrame::Ptr& key); + KeyFrameSnapshot(const Eigen::Isometry3d& pose, const pcl::PointCloud::ConstPtr& cloud); + + ~KeyFrameSnapshot(); + +public: + Eigen::Isometry3d pose; // pose estimated by graph optimization + pcl::PointCloud::ConstPtr cloud; // point cloud +}; + +} // namespace hdl_graph_slam + +#endif // KEYFRAME_HPP diff --git a/include/hdl_graph_slam/keyframe_updater.hpp b/include/hdl_graph_slam/keyframe_updater.hpp new file mode 100644 index 0000000..4f88e1e --- /dev/null +++ b/include/hdl_graph_slam/keyframe_updater.hpp @@ -0,0 +1,77 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef KEYFRAME_UPDATER_HPP +#define KEYFRAME_UPDATER_HPP + +#include +#include + +namespace hdl_graph_slam { + +/** + * @brief this class decides if a new frame should be registered to the pose graph as a keyframe + */ +class KeyframeUpdater { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + + /** + * @brief constructor + * @param pnh + */ + KeyframeUpdater(ros::NodeHandle& pnh) : is_first(true), prev_keypose(Eigen::Isometry3d::Identity()) { + keyframe_delta_trans = pnh.param("keyframe_delta_trans", 2.0); + keyframe_delta_angle = pnh.param("keyframe_delta_angle", 2.0); + + accum_distance = 0.0; + } + + /** + * @brief decide if a new frame should be registered to the graph + * @param pose pose of the frame + * @return if true, the frame should be registered + */ + bool update(const Eigen::Isometry3d& pose) { + // first frame is always registered to the graph + if(is_first) { + is_first = false; + prev_keypose = pose; + return true; + } + + // calculate the delta transformation from the previous keyframe + Eigen::Isometry3d delta = prev_keypose.inverse() * pose; + double dx = delta.translation().norm(); + double da = Eigen::AngleAxisd(delta.linear()).angle(); + + // too close to the previous frame + if(dx < keyframe_delta_trans && da < keyframe_delta_angle) { + return false; + } + + accum_distance += dx; + prev_keypose = pose; + return true; + } + + /** + * @brief the last keyframe's accumulated distance from the first keyframe + * @return accumulated distance + */ + double get_accum_distance() const { + return accum_distance; + } + +private: + // parameters + double keyframe_delta_trans; // + double keyframe_delta_angle; // + + bool is_first; + double accum_distance; + Eigen::Isometry3d prev_keypose; +}; + +} // namespace hdl_graph_slam + +#endif // KEYFRAME_UPDATOR_HPP diff --git a/include/hdl_graph_slam/loop_detector.hpp b/include/hdl_graph_slam/loop_detector.hpp new file mode 100644 index 0000000..3fd6b72 --- /dev/null +++ b/include/hdl_graph_slam/loop_detector.hpp @@ -0,0 +1,189 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef LOOP_DETECTOR_HPP +#define LOOP_DETECTOR_HPP + +#include +#include +#include +#include + +#include + +namespace hdl_graph_slam { + +struct Loop { +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + using Ptr = std::shared_ptr; + + Loop(const KeyFrame::Ptr& key1, const KeyFrame::Ptr& key2, const Eigen::Matrix4f& relpose) : key1(key1), key2(key2), relative_pose(relpose) {} + +public: + KeyFrame::Ptr key1; + KeyFrame::Ptr key2; + Eigen::Matrix4f relative_pose; +}; + +/** + * @brief this class finds loops by scam matching and adds them to the pose graph + */ +class LoopDetector { +public: + typedef pcl::PointXYZI PointT; + + /** + * @brief constructor + * @param pnh + */ + LoopDetector(ros::NodeHandle& pnh) { + distance_thresh = pnh.param("distance_thresh", 5.0); + accum_distance_thresh = pnh.param("accum_distance_thresh", 8.0); + distance_from_last_edge_thresh = pnh.param("min_edge_interval", 5.0); + + fitness_score_max_range = pnh.param("fitness_score_max_range", std::numeric_limits::max()); + fitness_score_thresh = pnh.param("fitness_score_thresh", 0.5); + + registration = select_registration_method(pnh); + last_edge_accum_distance = 0.0; + } + + /** + * @brief detect loops and add them to the pose graph + * @param keyframes keyframes + * @param new_keyframes newly registered keyframes + * @param graph_slam pose graph + */ + std::vector detect(const std::vector& keyframes, const std::deque& new_keyframes, hdl_graph_slam::GraphSLAM& graph_slam) { + std::vector detected_loops; + for(const auto& new_keyframe : new_keyframes) { + auto candidates = find_candidates(keyframes, new_keyframe); + auto loop = matching(candidates, new_keyframe, graph_slam); + if(loop) { + detected_loops.push_back(loop); + } + } + return detected_loops; + } + + double get_distance_thresh() const { + return distance_thresh; + } + +private: + /** + * @brief find loop candidates. A detected loop begins at one of #keyframes and ends at #new_keyframe + * @param keyframes candidate keyframes of loop start + * @param new_keyframe loop end keyframe + * @return loop candidates + */ + std::vector find_candidates(const std::vector& keyframes, const KeyFrame::Ptr& new_keyframe) const { + // too close to the last registered loop edge + std::cout<<"new_keyframe->accum_distance and distance "<accum_distance<<" "<accum_distance - last_edge_accum_distance < distance_from_last_edge_thresh) { + return std::vector(); + } + + std::vector candidates; + candidates.reserve(32); + + for(const auto& k : keyframes) { + // traveled distance between keyframes is too small + if(new_keyframe->accum_distance - k->accum_distance < accum_distance_thresh) { + continue; + } + + const auto& pos1 = k->node->estimate().translation(); + const auto& pos2 = new_keyframe->node->estimate().translation(); + + // estimated distance between keyframes is too small + double dist = (pos1.head<2>() - pos2.head<2>()).norm(); + if(dist > distance_thresh) { + continue; + } + + candidates.push_back(k); + } + + return candidates; + } + + /** + * @brief To validate a loop candidate this function applies a scan matching between keyframes consisting the loop. If they are matched well, the loop is added to the pose graph + * @param candidate_keyframes candidate keyframes of loop start + * @param new_keyframe loop end keyframe + * @param graph_slam graph slam + */ + Loop::Ptr matching(const std::vector& candidate_keyframes, const KeyFrame::Ptr& new_keyframe, hdl_graph_slam::GraphSLAM& graph_slam) { + if(candidate_keyframes.empty()) { + return nullptr; + } + + registration->setInputTarget(new_keyframe->cloud); + + double best_score = std::numeric_limits::max(); + KeyFrame::Ptr best_matched; + Eigen::Matrix4f relative_pose; + + std::cout << std::endl; + std::cout << "--- loop detection ---" << std::endl; + std::cout << "num_candidates: " << candidate_keyframes.size() << std::endl; + std::cout << "matching" << std::flush; + auto t1 = ros::Time::now(); + + pcl::PointCloud::Ptr aligned(new pcl::PointCloud()); + for(const auto& candidate : candidate_keyframes) { + registration->setInputSource(candidate->cloud); + Eigen::Isometry3d new_keyframe_estimate = new_keyframe->node->estimate(); + new_keyframe_estimate.linear() = Eigen::Quaterniond(new_keyframe_estimate.linear()).normalized().toRotationMatrix(); + Eigen::Isometry3d candidate_estimate = candidate->node->estimate(); + candidate_estimate.linear() = Eigen::Quaterniond(candidate_estimate.linear()).normalized().toRotationMatrix(); + Eigen::Matrix4f guess = (new_keyframe_estimate.inverse() * candidate_estimate).matrix().cast(); + guess(2, 3) = 0.0; + registration->align(*aligned, guess); + std::cout << "." << std::flush; + + double score = registration->getFitnessScore(fitness_score_max_range); + if(!registration->hasConverged() || score > best_score) { + continue; + } + + best_score = score; + best_matched = candidate; + relative_pose = registration->getFinalTransformation(); + } + + auto t2 = ros::Time::now(); + std::cout << " done" << std::endl; + std::cout << "best_score: " << boost::format("%.3f") % best_score << " time: " << boost::format("%.3f") % (t2 - t1).toSec() << "[sec]" << std::endl; + + if(best_score > fitness_score_thresh) { + std::cout << "loop not found..." << std::endl; + return nullptr; + } + + std::cout << "loop found!!" << std::endl; + std::cout << "relpose: " << relative_pose.block<3, 1>(0, 3) << " - " << Eigen::Quaternionf(relative_pose.block<3, 3>(0, 0)).coeffs().transpose() << std::endl; + + last_edge_accum_distance = best_matched->accum_distance; + + + return std::make_shared(new_keyframe, best_matched, relative_pose); + } + +private: + double distance_thresh; // estimated distance between keyframes consisting a loop must be less than this distance + double accum_distance_thresh; // traveled distance between ... + double distance_from_last_edge_thresh; // a new loop edge must far from the last one at least this distance + + double fitness_score_max_range; // maximum allowable distance between corresponding points + double fitness_score_thresh; // threshold for scan matching + + double last_edge_accum_distance; + + pcl::Registration::Ptr registration; +}; + +} // namespace hdl_graph_slam + +#endif // LOOP_DETECTOR_HPP diff --git a/include/hdl_graph_slam/map_cloud_generator.hpp b/include/hdl_graph_slam/map_cloud_generator.hpp new file mode 100644 index 0000000..b450f7c --- /dev/null +++ b/include/hdl_graph_slam/map_cloud_generator.hpp @@ -0,0 +1,34 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef MAP_CLOUD_GENERATOR_HPP +#define MAP_CLOUD_GENERATOR_HPP + +#include +#include +#include +#include + +namespace hdl_graph_slam { + +/** + * @brief this class generates a map point cloud from registered keyframes + */ +class MapCloudGenerator { +public: + using PointT = pcl::PointXYZI; + + MapCloudGenerator(); + ~MapCloudGenerator(); + + /** + * @brief generates a map point cloud + * @param keyframes snapshots of keyframes + * @param resolution resolution of generated map + * @return generated map point cloud + */ + pcl::PointCloud::Ptr generate(const std::vector& keyframes, double resolution) const; +}; + +} // namespace hdl_graph_slam + +#endif // MAP_POINTCLOUD_GENERATOR_HPP diff --git a/include/hdl_graph_slam/nmea_sentence_parser.hpp b/include/hdl_graph_slam/nmea_sentence_parser.hpp new file mode 100644 index 0000000..5133355 --- /dev/null +++ b/include/hdl_graph_slam/nmea_sentence_parser.hpp @@ -0,0 +1,108 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef NMEA_SENTENCE_PARSER_HPP +#define NMEA_SENTENCE_PARSER_HPP + +#include +#include +#include +#include +#include + +namespace hdl_graph_slam { + +struct GPRMC { +public: + GPRMC() { + status = 'V'; + } + + GPRMC(const std::vector& tokens) { + if(tokens[0] != "$GPRMC" || tokens.size() < 12) { + status = 'V'; + return; + } + + long time = std::stol(tokens[1]); + hour = time / 10000; + minute = (time % 10000) / 100; + second = time % 100; + + status = tokens[2][0]; + + latitude = degmin2deg(std::stod(tokens[3])); + latitude = tokens[4] == "N" ? latitude : -latitude; + + longitude = degmin2deg(std::stod(tokens[5])); + longitude = tokens[6] == "E" ? longitude : -longitude; + + speed_knots = std::stod(tokens[7]); + track_angle_degree = std::stod(tokens[8]); + + long date = std::stol(tokens[9]); + year = date % 100; + month = (date / 100) % 100; + day = (date / 10000) % 100; + + magnetic_variation = std::stod(tokens[10]); + magnetic_variation = tokens[11][0] == 'E' ? magnetic_variation : -magnetic_variation; + } + + double degmin2deg(double degmin) { + double d = std::floor(degmin / 100.0); + double m = (degmin - d * 100.0) / 60.0; + return d + m; + } + +public: + char status; // Status A=active or V=Void. + + int hour; // Fix taken at 12:35:19 UTC + int minute; + int second; + + double latitude; // + double longitude; + + double speed_knots; // Speed over the ground in knots + double track_angle_degree; // Track angle in degrees True + + int year; + int month; + int day; + + double magnetic_variation; +}; + +class NmeaSentenceParser { +public: + NmeaSentenceParser() {} + ~NmeaSentenceParser() {} + + GPRMC parse(const std::string& sentence) const { + int checksum_loc = sentence.find('*'); + if(checksum_loc == std::string::npos) { + return GPRMC(); + } + + int checksum = std::stoul(sentence.substr(checksum_loc + 1), nullptr, 16); + + std::string substr = sentence.substr(1, checksum_loc - 1); + int sum = std::accumulate(substr.begin(), substr.end(), static_cast(0), [=](unsigned char n, unsigned char c) { return n ^ c; }); + + if(checksum != (sum & 0xf)) { + std::cerr << "checksum doesn't match!!" << std::endl; + std::cerr << sentence << " " << sum << std::endl; + return GPRMC(); + } + + std::vector tokens; + boost::split(tokens, sentence, boost::is_any_of(",")); + + return GPRMC(tokens); + } +}; + +} // namespace hdl_graph_slam + +#endif // NMEA_SENTENCE_PARSER_HPP diff --git a/include/hdl_graph_slam/registrations.hpp b/include/hdl_graph_slam/registrations.hpp new file mode 100644 index 0000000..07af3bf --- /dev/null +++ b/include/hdl_graph_slam/registrations.hpp @@ -0,0 +1,21 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef HDL_GRAPH_SLAM_REGISTRATIONS_HPP +#define HDL_GRAPH_SLAM_REGISTRATIONS_HPP + +#include + +#include + +namespace hdl_graph_slam { + +/** + * @brief select a scan matching algorithm according to rosparams + * @param pnh + * @return selected scan matching + */ +pcl::Registration::Ptr select_registration_method(ros::NodeHandle& pnh); + +} // namespace hdl_graph_slam + +#endif // diff --git a/include/hdl_graph_slam/ros_time_hash.hpp b/include/hdl_graph_slam/ros_time_hash.hpp new file mode 100644 index 0000000..066085c --- /dev/null +++ b/include/hdl_graph_slam/ros_time_hash.hpp @@ -0,0 +1,24 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef ROS_TIME_HASH_HPP +#define ROS_TIME_HASH_HPP + +#include +#include + +#include + +/** + * @brief Hash calculation for ros::Time + */ +class RosTimeHash { +public: + size_t operator()(const ros::Time& val) const { + size_t seed = 0; + boost::hash_combine(seed, val.sec); + boost::hash_combine(seed, val.nsec); + return seed; + } +}; + +#endif // ROS_TIME_HASHER_HPP diff --git a/include/hdl_graph_slam/ros_utils.hpp b/include/hdl_graph_slam/ros_utils.hpp new file mode 100644 index 0000000..0a53617 --- /dev/null +++ b/include/hdl_graph_slam/ros_utils.hpp @@ -0,0 +1,93 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#ifndef ROS_UTILS_HPP +#define ROS_UTILS_HPP + +#include + +#include +#include +#include +#include + +namespace hdl_graph_slam { + +/** + * @brief convert Eigen::Matrix to geometry_msgs::TransformStamped + * @param stamp timestamp + * @param pose Eigen::Matrix to be converted + * @param frame_id tf frame_id + * @param child_frame_id tf child frame_id + * @return converted TransformStamped + */ +static geometry_msgs::TransformStamped matrix2transform(const ros::Time& stamp, const Eigen::Matrix4f& pose, const std::string& frame_id, const std::string& child_frame_id) { + Eigen::Quaternionf quat(pose.block<3, 3>(0, 0)); + quat.normalize(); + geometry_msgs::Quaternion odom_quat; + odom_quat.w = quat.w(); + odom_quat.x = quat.x(); + odom_quat.y = quat.y(); + odom_quat.z = quat.z(); + + geometry_msgs::TransformStamped odom_trans; + odom_trans.header.stamp = stamp; + odom_trans.header.frame_id = frame_id; + odom_trans.child_frame_id = child_frame_id; + + odom_trans.transform.translation.x = pose(0, 3); + odom_trans.transform.translation.y = pose(1, 3); + odom_trans.transform.translation.z = pose(2, 3); + odom_trans.transform.rotation = odom_quat; + + return odom_trans; +} + +static Eigen::Isometry3d pose2isometry(const geometry_msgs::Pose& pose) { + Eigen::Isometry3d mat = Eigen::Isometry3d::Identity(); + mat.translation() = Eigen::Vector3d(pose.position.x, pose.position.y, pose.position.z); + mat.linear() = Eigen::Quaterniond(pose.orientation.w, pose.orientation.x, pose.orientation.y, pose.orientation.z).toRotationMatrix(); + return mat; +} + +static Eigen::Isometry3d tf2isometry(const tf::StampedTransform& trans) { + Eigen::Isometry3d mat = Eigen::Isometry3d::Identity(); + mat.translation() = Eigen::Vector3d(trans.getOrigin().x(), trans.getOrigin().y(), trans.getOrigin().z()); + mat.linear() = Eigen::Quaterniond(trans.getRotation().w(), trans.getRotation().x(), trans.getRotation().y(), trans.getRotation().z()).toRotationMatrix(); + return mat; +} + +static geometry_msgs::Pose isometry2pose(const Eigen::Isometry3d& mat) { + Eigen::Quaterniond quat(mat.linear()); + Eigen::Vector3d trans = mat.translation(); + + geometry_msgs::Pose pose; + pose.position.x = trans.x(); + pose.position.y = trans.y(); + pose.position.z = trans.z(); + pose.orientation.w = quat.w(); + pose.orientation.x = quat.x(); + pose.orientation.y = quat.y(); + pose.orientation.z = quat.z(); + + return pose; +} + +static Eigen::Isometry3d odom2isometry(const nav_msgs::OdometryConstPtr& odom_msg) { + const auto& orientation = odom_msg->pose.pose.orientation; + const auto& position = odom_msg->pose.pose.position; + + Eigen::Quaterniond quat; + quat.w() = orientation.w; + quat.x() = orientation.x; + quat.y() = orientation.y; + quat.z() = orientation.z; + + Eigen::Isometry3d isometry = Eigen::Isometry3d::Identity(); + isometry.linear() = quat.toRotationMatrix(); + isometry.translation() = Eigen::Vector3d(position.x, position.y, position.z); + return isometry; +} + +} // namespace hdl_graph_slam + +#endif // ROS_UTILS_HPP diff --git a/include/hdl_graph_slam/utility.h b/include/hdl_graph_slam/utility.h new file mode 100644 index 0000000..f9d9d5f --- /dev/null +++ b/include/hdl_graph_slam/utility.h @@ -0,0 +1,174 @@ +#ifndef _UTILITY_LIDAR_ODOMETRY_H_ +#define _UTILITY_LIDAR_ODOMETRY_H_ + +#pragma once +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +typedef pcl::PointXYZI PointType; + +const string pointCloudTopic = "/kitti/velo/pointcloud"; +const string imuTopic = "/imu/data"; + + +// Save pcd +const string fileDirectory = "/home/wjk/3D Lidar_SLAM/data/"; + +// VLP-16 + +const int N_SCAN = 16; +const int Horizon_SCAN = 1800; +const float ang_res_x = 0.2; +const float ang_res_y = 2.0; +const float ang_bottom = 15.0+0.1; +const int groundScanInd = 7; + +// HDL-32E +/* + const int N_SCAN = 32; + const int Horizon_SCAN = 1800; + const float ang_res_x = 360.0/float(Horizon_SCAN); + const float ang_res_y = 41.33/float(N_SCAN-1); + const float ang_bottom = 30.67; + const int groundScanInd = 7;*/ + +// Ouster users may need to uncomment line 159 in imageProjection.cpp +// Usage of Ouster imu data is not fully supported yet, please just publish point cloud data +// Ouster OS1-16 +// const int N_SCAN = 16; +// const int Horizon_SCAN = 1024; +// const float ang_res_x = 360.0/float(Horizon_SCAN); +// const float ang_res_y = 33.2/float(N_SCAN-1); +// const float ang_bottom = 16.6+0.1; +// const int groundScanInd = 7; + +// Ouster OS1-64 +/* + const int N_SCAN = 64; + const int Horizon_SCAN = 1024; + const float ang_res_x = 360.0/float(Horizon_SCAN); + const float ang_res_y = 33.2/float(N_SCAN-1); + const float ang_bottom = 16.6+0.1; + const int groundScanInd = 15;*/ + +//Vel 64 +/* +const int N_SCAN = 64; +const int Horizon_SCAN = 1800; +const float ang_res_x = 0.2; +const float ang_res_y = 0.427; +const float ang_bottom = 24.9; +const int groundScanInd = 50;*/ + +//bool loopClosureEnableFlag; +const double mappingProcessInterval = 0.3; + +const float scanPeriod = 0.1; +const int systemDelay = 0; +const int imuQueLength = 200; + +const float sensorMountAngle = 0.0; +const float segmentTheta = 60.0/180.0*M_PI; // decrese this value may improve accuracy +const int segmentValidPointNum = 5; +const int segmentValidLineNum = 3; +const float segmentAlphaX = ang_res_x / 180.0 * M_PI; +const float segmentAlphaY = ang_res_y / 180.0 * M_PI; + + + const int edgeFeatureNum = 2; + const int surfFeatureNum = 4; + const int sectionsTotal = 6; + const float edgeThreshold = 0.1; + const float surfThreshold = 0.1; + const float nearestFeatureSearchSqDist = 25; + + + +// Mapping Params + const float surroundingKeyframeSearchRadius = 50.0; // key frame that is within n meters from current pose will be considerd for scan-to-map optimization (when loop closure disabled) + const int surroundingKeyframeSearchNum = 10; // submap size (when loop closure enabled) +// history key frames (history submap for loop closure) + const float historyKeyframeSearchRadius = 7.0; // key frame that is within n meters from current pose will be considerd for loop closure + const int historyKeyframeSearchNum = 25; // 2n+1 number of hostory key frames will be fused into a submap for loop closure + const float historyKeyframeFitnessScore = 0.3; // the smaller the better alignment + + const float globalMapVisualizationSearchRadius = 500.0; // key frames with in n meters will be visualized + + +struct smoothness_t{ + float value; + size_t ind; +}; + +struct by_value{ + bool operator()(smoothness_t const &left, smoothness_t const &right) { + return left.value < right.value; + } +}; + +/* + * A point cloud type that has 6D pose info ([x,y,z,roll,pitch,yaw] intensity is time stamp) + */ +struct PointXYZIRPYT +{ + PCL_ADD_POINT4D + PCL_ADD_INTENSITY; + float roll; + float pitch; + float yaw; + double time; + EIGEN_MAKE_ALIGNED_OPERATOR_NEW +} EIGEN_ALIGN16; + +POINT_CLOUD_REGISTER_POINT_STRUCT (PointXYZIRPYT, + (float, x, x) (float, y, y) + (float, z, z) (float, intensity, intensity) + (float, roll, roll) (float, pitch, pitch) (float, yaw, yaw) + (double, time, time) +) + +typedef PointXYZIRPYT PointTypePose; + +struct cloud_info{ + + int* startRingIndex; + int* endRingIndex; + + float startOrientation; + float endOrientation; + float orientationDiff; + + bool* segmentedCloudGroundFlag ; + int* segmentedCloudColInd; + float* segmentedCloudRange; +}; + +#endif diff --git a/launch/hdl_graph_slam.launch b/launch/hdl_graph_slam.launch new file mode 100644 index 0000000..b5e7868 --- /dev/null +++ b/launch/hdl_graph_slam.launch @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/hdl_graph_slam_400.launch b/launch/hdl_graph_slam_400.launch new file mode 100644 index 0000000..1b4b4f2 --- /dev/null +++ b/launch/hdl_graph_slam_400.launch @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/hdl_graph_slam_400_ours.launch b/launch/hdl_graph_slam_400_ours.launch new file mode 100644 index 0000000..27b8485 --- /dev/null +++ b/launch/hdl_graph_slam_400_ours.launch @@ -0,0 +1,165 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/hdl_graph_slam_501.launch b/launch/hdl_graph_slam_501.launch new file mode 100644 index 0000000..8f009c0 --- /dev/null +++ b/launch/hdl_graph_slam_501.launch @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/hdl_graph_slam_imu.launch b/launch/hdl_graph_slam_imu.launch new file mode 100644 index 0000000..b7c4119 --- /dev/null +++ b/launch/hdl_graph_slam_imu.launch @@ -0,0 +1,160 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/hdl_graph_slam_kitti.launch b/launch/hdl_graph_slam_kitti.launch new file mode 100644 index 0000000..b3b0f8d --- /dev/null +++ b/launch/hdl_graph_slam_kitti.launch @@ -0,0 +1,147 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/hdl_graph_slam_kitti_ours.launch b/launch/hdl_graph_slam_kitti_ours.launch new file mode 100644 index 0000000..ad0fef0 --- /dev/null +++ b/launch/hdl_graph_slam_kitti_ours.launch @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/hdl_graph_slam_mulran_ours.launch b/launch/hdl_graph_slam_mulran_ours.launch new file mode 100644 index 0000000..53b2bed --- /dev/null +++ b/launch/hdl_graph_slam_mulran_ours.launch @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/hdl_graph_slam_nclt_ours.launch b/launch/hdl_graph_slam_nclt_ours.launch new file mode 100644 index 0000000..1dba576 --- /dev/null +++ b/launch/hdl_graph_slam_nclt_ours.launch @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/hdl_graph_slam_stevens_ours.launch b/launch/hdl_graph_slam_stevens_ours.launch new file mode 100644 index 0000000..27b8485 --- /dev/null +++ b/launch/hdl_graph_slam_stevens_ours.launch @@ -0,0 +1,165 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/launch/msf_config.yaml b/launch/msf_config.yaml new file mode 100644 index 0000000..8206cc1 --- /dev/null +++ b/launch/msf_config.yaml @@ -0,0 +1,39 @@ +/msf_lidar_scan_matching/core/data_playback: true +/msf_lidar_scan_matching/core/fixed_bias: 0 + +#########IMU PARAMETERS####### +######## Asctec Firefly IMU +/msf_lidar_scan_matching/core/core_noise_acc: 0.083 +/msf_lidar_scan_matching/core/core_noise_accbias: 0.0083 +/msf_lidar_scan_matching/core/core_noise_gyr: 0.0013 +/msf_lidar_scan_matching/core/core_noise_gyrbias: 0.00013 + +/msf_lidar_scan_matching/pose_sensor/fixed_scale: 1 +/msf_lidar_scan_matching/pose_sensor/pose_noise_scale: 0.0 +/msf_lidar_scan_matching/pose_sensor/pose_noise_p_wv: 0.0 +/msf_lidar_scan_matching/pose_sensor/pose_noise_q_wv: 0.0 +/msf_lidar_scan_matching/pose_sensor/pose_noise_q_ic: 0.0 +/msf_lidar_scan_matching/pose_sensor/pose_noise_p_ic: 0.0 +/msf_lidar_scan_matching/pose_sensor/pose_delay: 0.02 +/msf_lidar_scan_matching/pose_sensor/pose_noise_meas_p: 0.05 +/msf_lidar_scan_matching/pose_sensor/pose_noise_meas_q: 0.05 + +/msf_lidar_scan_matching/pose_sensor/init/q_ic/w: 1.0 +/msf_lidar_scan_matching/pose_sensor/init/q_ic/x: 0.0 +/msf_lidar_scan_matching/pose_sensor/init/q_ic/y: 0.0 +/msf_lidar_scan_matching/pose_sensor/init/q_ic/z: 0.0 +/msf_lidar_scan_matching/pose_sensor/init/p_ic/x: 0.0 +/msf_lidar_scan_matching/pose_sensor/init/p_ic/y: 0.0 +/msf_lidar_scan_matching/pose_sensor/init/p_ic/z: 0.0 + +/msf_lidar_scan_matching/pose_sensor/pose_fixed_scale: true +/msf_lidar_scan_matching/pose_sensor/pose_fixed_p_ic: false +/msf_lidar_scan_matching/pose_sensor/pose_fixed_q_ic: false +/msf_lidar_scan_matching/pose_sensor/pose_fixed_p_wv: false +/msf_lidar_scan_matching/pose_sensor/pose_fixed_q_wv: false + +/msf_lidar_scan_matching/pose_sensor/pose_absolute_measurements: false +/msf_lidar_scan_matching/pose_sensor/pose_use_fixed_covariance: true +/msf_lidar_scan_matching/pose_sensor/pose_measurement_world_sensor: true # selects if sensor measures its position w.r.t. world (true, e.g. Vicon) or the position of the world coordinate system w.r.t. the sensor (false, e.g. ethzasl_ptam + +/msf_lidar_scan_matching/pose_sensor/pose_measurement_minimum_dt: 0.05 # Sets the minimum time in seconds between two pose measurements. diff --git a/msg/FloorCoeffs.msg b/msg/FloorCoeffs.msg new file mode 100644 index 0000000..514d53b --- /dev/null +++ b/msg/FloorCoeffs.msg @@ -0,0 +1,3 @@ +Header header + +float32[] coeffs diff --git a/msg/ScanMatchingStatus.msg b/msg/ScanMatchingStatus.msg new file mode 100644 index 0000000..58d8b52 --- /dev/null +++ b/msg/ScanMatchingStatus.msg @@ -0,0 +1,9 @@ +Header header + +bool has_converged +float32 matching_error +float32 inlier_fraction +geometry_msgs/Pose relative_pose + +std_msgs/String[] prediction_labels +geometry_msgs/Pose[] prediction_errors \ No newline at end of file diff --git a/nodelet_plugins.xml b/nodelet_plugins.xml new file mode 100644 index 0000000..f457ce5 --- /dev/null +++ b/nodelet_plugins.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/package.xml b/package.xml new file mode 100644 index 0000000..99203af --- /dev/null +++ b/package.xml @@ -0,0 +1,44 @@ + + + hdl_graph_slam + 0.0.0 + The hdl_graph_slam package + + koide + + BSD + + catkin + ndt_omp + fast_gicp + pcl_ros + roscpp + rospy + nodelet + geodesy + nmea_msgs + sensor_msgs + geometry_msgs + message_generation + libg2o + + ndt_omp + fast_gicp + pcl_ros + roscpp + rospy + geodesy + nodelet + nmea_msgs + sensor_msgs + message_generation + + python-scipy + python3-scipy + python-progressbar + python3-progressbar + + + + + diff --git a/setup.py b/setup.py new file mode 100755 index 0000000..8097795 --- /dev/null +++ b/setup.py @@ -0,0 +1,12 @@ +## ! DO NOT MANUALLY INVOKE THIS setup.py, USE CATKIN INSTEAD + +from setuptools import setup +from catkin_pkg.python_setup import generate_distutils_setup + +# fetch values from package.xml +setup_args = generate_distutils_setup( + packages=['hdl_graph_slam'], + package_dir={'': 'src'}) + +setup(**setup_args) + diff --git a/src/g2o/robust_kernel_io.cpp b/src/g2o/robust_kernel_io.cpp new file mode 100644 index 0000000..bd1c2e3 --- /dev/null +++ b/src/g2o/robust_kernel_io.cpp @@ -0,0 +1,156 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include + +#include +#include +#include +#include +#include +#include + +namespace g2o { + +std::string kernel_type(g2o::RobustKernel* kernel) { + if(dynamic_cast(kernel)) { + return "Huber"; + } + if(dynamic_cast(kernel)) { + return "Cauchy"; + } + if(dynamic_cast(kernel)) { + return "DCS"; + } + if(dynamic_cast(kernel)) { + return "Fair"; + } + if(dynamic_cast(kernel)) { + return "GemanMcClure"; + } + if(dynamic_cast(kernel)) { + return "PseudoHuber"; + } + if(dynamic_cast(kernel)) { + return "Saturated"; + } + if(dynamic_cast(kernel)) { + return "Tukey"; + } + if(dynamic_cast(kernel)) { + return "Welsch"; + } + return ""; +} + +bool save_robust_kernels(const std::string& filename, g2o::SparseOptimizer* graph) { + std::ofstream ofs(filename); + if(!ofs) { + std::cerr << "failed to open output stream!!" << std::endl; + return false; + } + + for(const auto& edge_ : graph->edges()) { + g2o::OptimizableGraph::Edge* edge = static_cast(edge_); + g2o::RobustKernel* kernel = edge->robustKernel(); + if(!kernel) { + continue; + } + + std::string type = kernel_type(kernel); + if(type.empty()) { + std::cerr << "unknown kernel type!!" << std::endl; + continue; + } + + ofs << edge->vertices().size() << " "; + for(size_t i = 0; i < edge->vertices().size(); i++) { + ofs << edge->vertices()[i]->id() << " "; + } + ofs << type << " " << kernel->delta() << std::endl; + } + + return true; +} + +class KernelData { +public: + KernelData(const std::string& line) { + std::stringstream sst(line); + size_t num_vertices; + sst >> num_vertices; + + vertex_indices.resize(num_vertices); + for(size_t i = 0; i < num_vertices; i++) { + sst >> vertex_indices[i]; + } + + sst >> type >> delta; + } + + bool match(g2o::OptimizableGraph::Edge* edge) const { + if(edge->vertices().size() != vertex_indices.size()) { + return false; + } + + for(size_t i = 0; i < edge->vertices().size(); i++) { + if(edge->vertices()[i]->id() != vertex_indices[i]) { + return false; + } + } + + return true; + } + + g2o::RobustKernel* create() const { + auto factory = g2o::RobustKernelFactory::instance(); + g2o::RobustKernel* kernel = factory->construct(type); + kernel->setDelta(delta); + + return kernel; + } + +public: + std::vector vertex_indices; + std::string type; + double delta; +}; + +bool load_robust_kernels(const std::string& filename, g2o::SparseOptimizer* graph) { + std::ifstream ifs(filename); + if(!ifs) { + std::cerr << "warning: failed to open input stream!!" << std::endl; + return true; + } + + std::vector kernels; + + while(!ifs.eof()) { + std::string line; + std::getline(ifs, line); + if(line.empty()) { + continue; + } + + kernels.push_back(KernelData(line)); + } + std::cout << "kernels: " << kernels.size() << std::endl; + + for(auto& edge_ : graph->edges()) { + g2o::OptimizableGraph::Edge* edge = static_cast(edge_); + + for(auto itr = kernels.begin(); itr != kernels.end(); itr++) { + if(itr->match(edge)) { + edge->setRobustKernel(itr->create()); + kernels.erase(itr); + break; + } + } + } + + if(kernels.size() != 0) { + std::cerr << "warning: there is non-associated kernels!!" << std::endl; + } + return true; +} + +} // namespace g2o diff --git a/src/hdl_graph_slam/__init__.py b/src/hdl_graph_slam/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/hdl_graph_slam/bag_player.py b/src/hdl_graph_slam/bag_player.py new file mode 100755 index 0000000..d161627 --- /dev/null +++ b/src/hdl_graph_slam/bag_player.py @@ -0,0 +1,223 @@ +#!/usr/bin/python +# SPDX-License-Identifier: BSD-2-Clause +import sys +import yaml +import time +import curses +import argparse + +try: + from StringIO import StringIO ## for Python 2 +except ImportError: + from io import StringIO ## for Python 3 + + +import rospy +import rosbag +import roslib +from std_msgs.msg import * +from sensor_msgs.msg import * +from rosgraph_msgs.msg import * +from progressbar import ProgressBar + + +class BagPlayer: + def __init__(self, bagfile, start, duration): + print('opening...',) + self.bag = rosbag.Bag(bagfile, 'r') + print('done') + + self.message_generator = self.bag.read_messages() + + info_dict = yaml.safe_load(self.bag._get_yaml_info()) + self.duration = float(info_dict['duration']) + self.endtime = float(info_dict['end']) + + self.progress = ProgressBar(0, self.duration, term_width=80, fd=StringIO()) + + self.publishers = {} + for con in self.bag._get_connections(): + msg_class = roslib.message.get_message_class(con.datatype) + self.publishers[con.topic] = rospy.Publisher(con.topic, msg_class, queue_size=256) + self.clock_pub = rospy.Publisher('/clock', Clock, queue_size=256) + + self.init_time = None + self.start = start + self.duration = duration + self.fail_count = 0 + self.time_subs = {} + self.target_times = {} + self.latest_stamps = {} + + self.play() + + def update_time_subs(self): + for topic_name, msg_type in rospy.get_published_topics(): + if 'read_until' in topic_name and 'std_msgs/Header' in msg_type: + if topic_name not in self.time_subs: + print('ADD', topic_name) + self.time_subs[topic_name] = rospy.Subscriber(topic_name, Header, self.time_callback, topic_name) + + def time_callback(self, header_msg, topic_name): + if header_msg.frame_id not in self.target_times: + self.target_times[header_msg.frame_id] = {} + + if topic_name not in self.target_times[header_msg.frame_id] or self.target_times[header_msg.frame_id][topic_name] < header_msg.stamp: + self.target_times[header_msg.frame_id][topic_name] = header_msg.stamp + + def play_realtime(self, duration): + topic, msg, stamp = next(self.message_generator) + stamp_wall = time.time() + + start_stamp = stamp + start_stamp_wall = stamp_wall + + self.start_stamp = start_stamp + + while not rospy.is_shutdown() and (stamp - start_stamp) < duration: + stamp_wall = time.time() + elapsed_stamp = stamp - start_stamp + if (stamp_wall - start_stamp_wall) < (elapsed_stamp.secs + elapsed_stamp.nsecs * 1e-9): + time.sleep(1e-6) + self.update_time_subs() + continue + + clock_msg = Clock() + clock_msg.clock = stamp + + if self.init_time is None: + self.init_time = stamp + + if self.start and (stamp - self.init_time) < rospy.Duration(float(self.start)): + start_stamp = stamp + else: + self.clock_pub.publish(clock_msg) + self.publishers[topic].publish(msg) + + topic, msg, stamp = next(self.message_generator) + + return topic, msg, stamp + + def print_progress(self, stamp): + self.stdscr.clear() + self.stdscr.addstr(0, 0, 'topic') + self.stdscr.addstr(0, 25, 'time') + + line = 1 + for target in self.target_times: + if target not in self.publishers: + continue + + for sub_name in self.target_times[target]: + target_time = self.target_times[target][sub_name] + self.stdscr.addstr(line, 0, sub_name[:-11]) + self.stdscr.addstr(line, 25, '%.6f' % (target_time.secs + target_time.nsecs * 1e-9)) + + residual = target_time - self.latest_stamps[target].stamp + + color = 1 if residual.to_sec() > 0.0 else 2 + self.stdscr.addstr(line, 50, '%.5f' % residual.to_sec(), curses.color_pair(color)) + line += 1 + + if not hasattr(self, 'prev_stamp'): + self.prev_stamp = stamp + self.prev_stamp_wall = time.time() + self.processing_speed = 1.0 + elif time.time() - self.prev_stamp_wall > 1.0: + sim_duration = (stamp - self.prev_stamp).to_sec() + wall_duration = time.time() - self.prev_stamp_wall + self.processing_speed = sim_duration / wall_duration + + self.stdscr.addstr(line, 0, 'current_stamp') + self.stdscr.addstr(line, 25, '%.6f' % stamp.to_sec()) + self.stdscr.addstr(line, 50, '(x%.2f)' % self.processing_speed) + + elapsed = (stamp - self.start_stamp).to_sec() + self.progress.fd = StringIO() + try: + self.progress.update(elapsed) + except: + # nothing to do + pass + self.stdscr.addstr(line + 1, 0, '----------') + self.stdscr.addstr(line + 2, 0, self.progress.fd.getvalue()) + + self.stdscr.refresh() + + def check_stamp(self, topic, msg): + if topic not in self.target_times: + return True + + if self.fail_count > 10: + self.fail_count = 0 + return True + + target_time_map = self.target_times[topic] + for sub_name in target_time_map: + self.latest_stamps[topic] = msg.header + if msg.header.stamp > target_time_map[sub_name]: + self.fail_count += 1 + return False + + self.fail_count = 0 + return True + + def play(self): + print('play realtime for 3.0[sec]') + topic, msg, stamp = self.play_realtime(rospy.Duration(5.0)) + self.update_time_subs() + + print('play as fast as possible') + self.stdscr = curses.initscr() + curses.start_color() + curses.noecho() + + curses.init_pair(1, curses.COLOR_BLUE, curses.COLOR_WHITE) + curses.init_pair(2, curses.COLOR_RED, curses.COLOR_WHITE) + + try: + while not rospy.is_shutdown(): + if not self.check_stamp(topic, msg): + self.update_time_subs() + self.print_progress(stamp) + time.sleep(0.1) + continue + + clock_msg = Clock() + clock_msg.clock = stamp + + if self.duration: + if (stamp - self.init_time) > rospy.Duration(float(self.duration)): + break + + self.clock_pub.publish(clock_msg) + self.publishers[topic].publish(msg) + topic, msg, stamp = next(self.message_generator) + except: + print(sys.exc_info()[0]) + clock_msg = Clock() + clock_msg.clock = stamp + rospy.Duration(30.0) + self.clock_pub.publish(clock_msg) + time.sleep(0.5) + + curses.echo() + curses.endwin() + + +def main(): + myargv = rospy.myargv(sys.argv) + parser = argparse.ArgumentParser() + parser.add_argument('input_bag', help='bag file to be played') + parser.add_argument('-s', '--start', help='start sec seconds into the bag') + parser.add_argument('-u', '--duration', help='play only sec seconds into the bag') + args = parser.parse_args(myargv[1:]) + + if len(sys.argv) < 2: + print('usage bag_player src_bagname') + return + + rospy.init_node('bag_player') + BagPlayer(args.input_bag, args.start, args.duration) + +if __name__ == '__main__': + main() diff --git a/src/hdl_graph_slam/featureAssociation.cpp b/src/hdl_graph_slam/featureAssociation.cpp new file mode 100644 index 0000000..0e25c39 --- /dev/null +++ b/src/hdl_graph_slam/featureAssociation.cpp @@ -0,0 +1,693 @@ +#include "hdl_graph_slam/featureAssociation.h" +#include "hdl_graph_slam/imageProjection.h" +#include +#include + +using namespace hdl_graph_slam; + +FeatureAssociation::FeatureAssociation() { + + initializationValue(); +} + +void FeatureAssociation::initializationValue() { + max_prob_ = 0.0; + avg_error_last_ = 0.2; + cloudSmoothness.resize(N_SCAN * Horizon_SCAN); + + downSizeFilter.setLeafSize(0.2, 0.2, 0.2); + + segmentedCloud.reset(new pcl::PointCloud()); + outlierCloud.reset(new pcl::PointCloud()); + + Source_filtered_.reset(new pcl::PointCloud()); + Target_filtered_.reset(new pcl::PointCloud()); + + LinearCloud_.makeShared().reset(new pcl::PointCloud()); + Less_LinearCloud_.makeShared().reset(new pcl::PointCloud()); + LinearCloud_Ground_.makeShared().reset(new pcl::PointCloud()); + LinearCloud_Less_Ground_.makeShared().reset(new pcl::PointCloud()); + + Less_LinearCloud_last_.makeShared().reset(new pcl::PointCloud()); + LinearCloud_Less_Ground_Last.makeShared().reset(new pcl::PointCloud()); + + cornerPointsSharp.reset(new pcl::PointCloud()); + cornerPointsLessSharp.reset(new pcl::PointCloud()); + surfPointsFlat.reset(new pcl::PointCloud()); + surfPointsLessFlat.reset(new pcl::PointCloud()); + + surfPointsLessFlatScan.reset(new pcl::PointCloud()); + surfPointsLessFlatScanDS.reset(new pcl::PointCloud()); + + timeScanCur = 0; + timeNewSegmentedCloud = 0; + timeNewSegmentedCloudInfo = 0; + timeNewOutlierCloud = 0; + + newSegmentedCloud = false; + newSegmentedCloudInfo = false; + newOutlierCloud = false; + + systemInitCount = 0; + systemInited = false; + + for (int i = 0; i < 6; ++i) { + transformCur[i] = 0; + transformSum[i] = 0; + } + + systemInitedLM = false; + + laserCloudCornerLast.reset(new pcl::PointCloud()); + laserCloudSurfLast.reset(new pcl::PointCloud()); + laserCloudOri.reset(new pcl::PointCloud()); + coeffSel.reset(new pcl::PointCloud()); + target_cloud_.reset(new pcl::PointCloud()); + + kdtreeLineCornerLast.reset(new pcl::KdTreeFLANN()); + kdtreeCornerLast.reset(new pcl::KdTreeFLANN()); + kdtreeSurfLast.reset(new pcl::KdTreeFLANN()); + kdtreeGroundCornerLast.reset(new pcl::KdTreeFLANN()); + + isDegenerate = false; + matP = cv::Mat(6, 6, CV_32F, cv::Scalar::all(0)); + + frameCount = 1; + + smooth_dist.setZero(); + smooth_dist_no_.setZero(); +} + +void FeatureAssociation::setSourceCloud(const pcl::PointCloud::ConstPtr cloud) +{ + pcl::copyPointCloud(*cloud,LaserCloud_); + + ImageProjection IP_source; + IP_source.Thre_ = Thre_; + IP_source.Time_ = Time_; + clock_t s1 = clock(); + IP_source.cloudHandler(LaserCloud_.makeShared()); + clock_t e1 = clock(); + pcl::copyPointCloud(*IP_source.LinearPoints_,LinearCloud_); + pcl::copyPointCloud(*IP_source.Less_LinearPoints_,Less_LinearCloud_); + pcl::copyPointCloud(*IP_source.LinearPoints_Less_Ground,LinearCloud_Less_Ground_); + pcl::copyPointCloud(*IP_source.LinearPoints_Ground,LinearCloud_Ground_); + clock_t e2 = clock(); + + cout<<"setSourceCloud "<<(double)(e2-e1)/CLOCKS_PER_SEC<<" "<<(double)(e1-s1)/CLOCKS_PER_SEC<::ConstPtr cloud) +{ + + pcl::copyPointCloud(*cloud,LaserCloud_last_); + /* + ImageProjection IP_target; + IP_target.Thre_ = Thre_; + IP_target.Time_ = Time_; + IP_target.cloudHandler(LaserCloud_last_.makeShared());*/ + pcl::copyPointCloud(LinearCloud_,LinearCloud_last_); + pcl::copyPointCloud(Less_LinearCloud_,Less_LinearCloud_last_); + pcl::copyPointCloud(LinearCloud_Less_Ground_,LinearCloud_Less_Ground_Last); + pcl::copyPointCloud(LinearCloud_Ground_,LinearCloud_Ground_Last); + cout<<"setTargetCloud"<::ConstPtr cloud) +{ + + pcl::copyPointCloud(*cloud,LaserCloud_last_); + + ImageProjection IP_target; + IP_target.Thre_ = Thre_; + IP_target.Time_ = Time_; + IP_target.cloudHandler(LaserCloud_last_.makeShared()); + pcl::copyPointCloud(*IP_target.LinearPoints_,LinearCloud_last_); + pcl::copyPointCloud(*IP_target.Less_LinearPoints_,Less_LinearCloud_last_); + pcl::copyPointCloud(*IP_target.LinearPoints_Less_Ground,LinearCloud_Less_Ground_Last); + pcl::copyPointCloud(*IP_target.LinearPoints_Ground,LinearCloud_Ground_Last); + cout<<"setTargetCloud"<intensity - int(pi->intensity)); + + float rx = s * transformCur[0]; + float ry = s * transformCur[1]; + float rz = s * transformCur[2]; + float tx = s * transformCur[3]; + float ty = s * transformCur[4]; + float tz = s * transformCur[5]; + + float x1 = cos(rz) * (pi->x) - sin(rz) * (pi->y); + float y1 = sin(rz) * (pi->x) + cos(rz) * (pi->y); + float z1 = (pi->z); + + float x2 = cos(ry) * x1 + sin(ry) * z1; + float y2 = y1; + float z2 = -sin(ry) * x1 + cos(ry) * z1; + + po->x = x2 + tx; + po->y = cos(rx) * y2 - sin(rx) * z2 + ty; + po->z = sin(rx) * y2 + cos(rx) * z2 + tz; + po->intensity = pi->intensity; +} + +double FeatureAssociation::rad2deg(double radians) { + return radians * 180.0 / M_PI; +} + +void FeatureAssociation::findCorrespondingCornerFeatures_X(int iterCount) { + smooth_dist.setOnes(); + smooth_dist_no_.setZero(); + total_cost_ = 0; + //std::cout<<"linear points: "<setInputCloud(Less_LinearCloud_last_.makeShared()); + int segmentedCloudNum = LinearCloud_.points.size(); + //std::cout<<"linear points: "<setInputCloud(laserCloudCornerLast); + PointType p0 = pointSel; + //p0.z = 0; + kdtreeLineCornerLast->nearestKSearch(p0, 1, pointSearchInd, pointSearchSqDis); + int closestPointInd = -1, minPointInd2 = -1; + // the two indexs + if (pointSearchSqDis[0] < nearestFeatureSearchSqDist) { + closestPointInd = pointSearchInd[0]; + } + + pointSearchCornerInd1[i] = closestPointInd; + } + //compute the point-to-line distance according to the cosine theory and the deviretion + if (pointSearchCornerInd1[i] >= 0) { + + tripod1 = Less_LinearCloud_last_.points[pointSearchCornerInd1[i]]; + //tripod2 = Less_LinearCloud_last_D_.points[pointSearchCornerInd1[i]]; + + float x0 = pointSel.x; + float y0 = pointSel.y; + float z0 = pointSel.z; + float x1 = tripod1.x; + float y1 = tripod1.y; + float z1 = tripod1.z; + float x2 = tripod1.x; + float y2 = tripod1.y; + float z2 = tripod1.z + 0.5; + + float m11 = ((x0 - x1) * (y0 - y2) - (x0 - x2) * (y0 - y1)); + float m22 = ((x0 - x1) * (z0 - z2) - (x0 - x2) * (z0 - z1)); + float m33 = ((y0 - y1) * (z0 - z2) - (y0 - y2) * (z0 - z1)); + + float a012 = sqrt(m11 * m11 + m22 * m22 + m33 * m33); + + float l12 = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2)); + + float la = ((y1 - y2) * m11 + (z1 - z2) * m22) / a012 / l12; + + float lb = -((x1 - x2) * m11 - (z1 - z2) * m33) / a012 / l12; + + float lc = -((x1 - x2) * m22 + (y1 - y2) * m33) / a012 / l12; + + float ld2 = a012 / l12; + + float s = 1; + if (iterCount >= 5) { + s = 1 - 1.8 * fabs(ld2); + } + + + error_sum += abs(ld2); + + if (s > 0.1) { + coeff.x = s * la; + coeff.y = s * lb; + coeff.z = s * lc; + coeff.intensity = s * ld2; + LinearCloud_.points[i].intensity = ld2; + laserCloudOri->push_back(LinearCloud_.points[i]); + coeffSel->push_back(coeff); + // + target_cloud_->push_back(tripod1); + } + } + } +} + + + + +void FeatureAssociation::findCorrespondingSurfFeatures(int iterCount) { + //cout<<"LinearCloud_Less_Ground_Last "<setInputCloud(LinearCloud_Less_Ground_Last.makeShared()); + total_cost_ = 0; + int surfPointsFlatNum = LinearCloud_Ground_.points.size(); + + for (int i = 0; i < surfPointsFlatNum; i++) { + //cout<<"begin transformto start"<nearestKSearch(pointSel, 1, pointSearchInd, pointSearchSqDis); + int closestPointInd = -1, minPointInd2 = -1, minPointInd3 = -1; + + if (pointSearchSqDis[0] < nearestFeatureSearchSqDist) { + closestPointInd = pointSearchInd[0]; + int closestPointScan = int(LinearCloud_Less_Ground_Last.points[closestPointInd].intensity); + + float pointSqDis, minPointSqDis2 = nearestFeatureSearchSqDist, minPointSqDis3 = nearestFeatureSearchSqDist; + for (int j = closestPointInd + 1; j < surfPointsFlatNum; j++) { + if (int(LinearCloud_Less_Ground_Last.points[j].intensity) > closestPointScan + 2.5) { + break; + } + + pointSqDis = (LinearCloud_Less_Ground_Last.points[j].x - pointSel.x) * + (LinearCloud_Less_Ground_Last.points[j].x - pointSel.x) + + (LinearCloud_Less_Ground_Last.points[j].y - pointSel.y) * + (LinearCloud_Less_Ground_Last.points[j].y - pointSel.y) + + (LinearCloud_Less_Ground_Last.points[j].z - pointSel.z) * + (LinearCloud_Less_Ground_Last.points[j].z - pointSel.z); + + if (int(LinearCloud_Less_Ground_Last.points[j].intensity) <= closestPointScan) { + if (pointSqDis < minPointSqDis2) { + minPointSqDis2 = pointSqDis; + minPointInd2 = j; + } + } else { + if (pointSqDis < minPointSqDis3) { + minPointSqDis3 = pointSqDis; + minPointInd3 = j; + } + } + } + for (int j = closestPointInd - 1; j >= 0; j--) { + if (int(LinearCloud_Less_Ground_Last.points[j].intensity) < closestPointScan - 2.5) { + break; + } + + pointSqDis = (LinearCloud_Less_Ground_Last.points[j].x - pointSel.x) * + (LinearCloud_Less_Ground_Last.points[j].x - pointSel.x) + + (LinearCloud_Less_Ground_Last.points[j].y - pointSel.y) * + (LinearCloud_Less_Ground_Last.points[j].y - pointSel.y) + + (LinearCloud_Less_Ground_Last.points[j].z - pointSel.z) * + (LinearCloud_Less_Ground_Last.points[j].z - pointSel.z); + + if (int(LinearCloud_Less_Ground_Last.points[j].intensity) >= closestPointScan) { + if (pointSqDis < minPointSqDis2) { + minPointSqDis2 = pointSqDis; + minPointInd2 = j; + } + } else { + if (pointSqDis < minPointSqDis3) { + minPointSqDis3 = pointSqDis; + minPointInd3 = j; + } + } + } + } + + pointSearchSurfInd1[i] = closestPointInd; + pointSearchSurfInd2[i] = minPointInd2; + pointSearchSurfInd3[i] = minPointInd3; + } + + if (pointSearchSurfInd2[i] >= 0 && pointSearchSurfInd3[i] >= 0) { + + tripod1 = LinearCloud_Less_Ground_Last.points[pointSearchSurfInd1[i]]; + tripod2 = LinearCloud_Less_Ground_Last.points[pointSearchSurfInd2[i]]; + tripod3 = LinearCloud_Less_Ground_Last.points[pointSearchSurfInd3[i]]; + + tripod2.z = tripod1.z; + tripod3.z = tripod1.z; + + float pa = (tripod2.y - tripod1.y) * (tripod3.z - tripod1.z) + - (tripod3.y - tripod1.y) * (tripod2.z - tripod1.z); + float pb = (tripod2.z - tripod1.z) * (tripod3.x - tripod1.x) + - (tripod3.z - tripod1.z) * (tripod2.x - tripod1.x); + float pc = (tripod2.x - tripod1.x) * (tripod3.y - tripod1.y) + - (tripod3.x - tripod1.x) * (tripod2.y - tripod1.y); + float pd = -(pa * tripod1.x + pb * tripod1.y + pc * tripod1.z); + + float ps = sqrt(pa * pa + pb * pb + pc * pc); + + pa /= ps; + pb /= ps; + pc /= ps; + pd /= ps; + + float pd2 = pa * pointSel.x + pb * pointSel.y + pc * pointSel.z + pd; + + float s = 1; + if (iterCount >= 5) { + s = 1 - 1.8 * fabs(pd2) / sqrt(sqrt(pointSel.x * pointSel.x + + pointSel.y * pointSel.y + pointSel.z * pointSel.z)); + } + + if (s > 0.1 && abs(pd2)>0) { + coeff.x = s * pa; + coeff.y = s * pb; + coeff.z = s * pc; + coeff.intensity = s * pd2; + total_cost_ = total_cost_ + pd2; + laserCloudOri->push_back(LinearCloud_Ground_.points[i]); + //tripod1.x = pointSel.x; + //tripod1.y = pointSel.y; + target_cloud_->push_back(tripod1); + coeffSel->push_back(coeff); + } + } + } +} + +bool FeatureAssociation::calculateTransformationSurf(int iterCount) { + + int pointSelNum = laserCloudOri->points.size(); + //cout<<"Ground pointSelNum: "<points[i]; + coeff = coeffSel->points[i]; + + float arx = ((crx*crz*sry-sry*srz)*pointOri.x + (-crz*srx-crx*sry*srz)*pointOri.y - crx*cry*pointOri.z)*coeff.y + + ((crx*srz+crz*srx*sry)*pointOri.x + (crx*crz-srx*sry*srz)*pointOri.y + (-cry*srx)*pointOri.z)*coeff.z; + + float ary = ((-crz*sry)*pointOri.x + sry*srz*pointOri.y + cry*pointOri.z)*coeff.x + + (cry*crz*srx*pointOri.x - cry*srx*srz*pointOri.y + srx*sry*pointOri.z)*coeff.y + + (-crx*cry*crz*pointOri.x + crx*cry*srz*pointOri.y - crx*sry*pointOri.z)*coeff.z; + + float atz = coeff.z; + + float d2 = coeff.intensity; + + matA.at(i, 0) = arx; + matA.at(i, 1) = ary; + matA.at(i, 2) = atz; + matB.at(i, 0) = -0.05 * d2; + } + + cv::transpose(matA, matAt); + matAtA = matAt * matA; + matAtB = matAt * matB; + //cout<<"begin solve"<= 0; i--) { + if (matE.at(0, i) < eignThre[i]) { + for (int j = 0; j < 3; j++) { + matV2.at(i, j) = 0; + } + isDegenerate = true; + } else { + break; + } + } + matP = matV.inv() * matV2; + } + + if (isDegenerate) { + cv::Mat matX2(3, 1, CV_32F, cv::Scalar::all(0)); + matX.copyTo(matX2); + matX = matP * matX2; + } + + transformCur[0] += matX.at(0, 0); + transformCur[1] += matX.at(1, 0); + transformCur[5] += matX.at(2, 0); + + for (int i = 0; i < 6; i++) { + if (isnan(transformCur[i])) + transformCur[i] = 0; + } + + float deltaR = sqrt( + pow(rad2deg(matX.at(0, 0)), 2) + + pow(rad2deg(matX.at(1, 0)), 2)); + float deltaT = sqrt( + pow(matX.at(2, 0) * 100, 2)); + + if (deltaR < 0.1 && deltaT < 0.1) { + return false; + } + return true; +} + +bool FeatureAssociation::calculateTransformationCorner(int iterCount) { + int pointSelNum = laserCloudOri->points.size(); + //cout<<"Corner pointSelNum: "<points[i]; + coeff = coeffSel->points[i]; + + float arz = (-cry*srz * pointOri.x - cry*crz * pointOri.y ) * coeff.x + + ((crx*crz-srx*sry*srz) * pointOri.x + (-crx*srz-crz*srx*sry) * pointOri.y ) * coeff.y + + ((crz*srx+crx*sry*srz)*pointOri.x + (crx*crz*sry-srx*srz)*pointOri.y)* coeff.z; + + float atx = coeff.x; + + float aty = coeff.y; + + float d2 = coeff.intensity; + + matA.at(i, 0) = arz; + matA.at(i, 1) = atx; + matA.at(i, 2) = aty; + matB.at(i, 0) = -0.1 * d2; + } + + cv::transpose(matA, matAt); + matAtA = matAt * matA; + matAtB = matAt * matB; + + cv::solve(matAtA, matAtB, matX, cv::DECOMP_QR); + + //sleep(1000); + if (iterCount == 0) { + cv::Mat matE(1, 3, CV_32F, cv::Scalar::all(0)); + cv::Mat matV(3, 3, CV_32F, cv::Scalar::all(0)); + cv::Mat matV2(3, 3, CV_32F, cv::Scalar::all(0)); + + cv::eigen(matAtA, matE, matV); + matV.copyTo(matV2); + + isDegenerate = false; + float eignThre[3] = {10, 10, 10}; + for (int i = 2; i >= 0; i--) { + if (matE.at(0, i) < eignThre[i]) { + for (int j = 0; j < 3; j++) { + matV2.at(i, j) = 0; + } + isDegenerate = true; + } else { + break; + } + } + matP = matV.inv() * matV2; + } + + if (isDegenerate) { + cv::Mat matX2(3, 1, CV_32F, cv::Scalar::all(0)); + matX.copyTo(matX2); + matX = matP * matX2; + } + transformCur[2] += matX.at(0, 0); + transformCur[3] += matX.at(1, 0); + transformCur[4] += matX.at(2, 0); + + for (int i = 0; i < 6; i++) { + if (isnan(transformCur[i])) { + transformCur[i] = 0; + } + + } + + float deltaR = sqrt( + pow(rad2deg(matX.at(0, 0)), 2)); + float deltaT = sqrt( + pow(matX.at(1, 0) * 100, 2) + + pow(matX.at(2, 0) * 100, 2)); + + if (deltaR < 0.1 && deltaT < 0.1) { + return false; + } + return true; +} + + +void FeatureAssociation::updateTransformation() { + iterCount1_ = 0; + iterCount2_ = 0; + + + for (int iterCount1 = 0; iterCount1 < 25; iterCount1++) { + laserCloudOri->clear(); + coeffSel->clear(); + + findCorrespondingSurfFeatures(iterCount1); + if (laserCloudOri->points.size() < 10) + continue; + if (calculateTransformationSurf(iterCount1) == false) + break; + iterCount1_ = iterCount1 + 1; + } + + for (int iterCount2 = 0; iterCount2 < 25; iterCount2++) { + laserCloudOri->clear(); + + coeffSel->clear(); + + findCorrespondingCornerFeatures_X(iterCount2); + + if (laserCloudOri->points.size() < 10) + continue; + + if (calculateTransformationCorner(iterCount2) == false) { + break; + } + } +} + +void FeatureAssociation::toEulerAngle( Eigen::Quaternionf& q, float& roll, float& pitch, float& yaw) +{ +// roll (x-axis rotation) + double sinr_cosp = +2.0 * (q.w() * q.x() + q.y() * q.z()); + double cosr_cosp = +1.0 - 2.0 * (q.x() * q.x() + q.y() * q.y()); + roll = atan2(sinr_cosp, cosr_cosp); + +// pitch (y-axis rotation) + double sinp = +2.0 * (q.w() * q.y() - q.z() * q.x()); + if (fabs(sinp) >= 1) + pitch = copysign(M_PI / 2, sinp); // use 90 degrees if out of range + else + pitch = asin(sinp); + +// yaw (z-axis rotation) + double siny_cosp = +2.0 * (q.w() * q.z() + q.x() * q.y()); + double cosy_cosp = +1.0 - 2.0 * (q.y() * q.y() + q.z() * q.z()); + yaw = atan2(siny_cosp, cosy_cosp); +} + +void FeatureAssociation::Matrix4dToRpyxyz(const Eigen::Matrix4f &Trans,float nums[]){ + + Eigen::Quaternionf quaternion(Trans.block<3,3>(0,0)); + float roll, pitch, yaw; + toEulerAngle(quaternion, roll, pitch, yaw); + + Eigen::Vector3f rpy_transformSum=Trans.block<3,3>(0,0).eulerAngles(2,1,0); + nums[0]= roll;//rpy_transformSum(0); + nums[1]= pitch;//rpy_transformSum(1);; + nums[2]= yaw;//rpy_transformSum(2);//; + nums[3]=Trans(0,3); + nums[4]=Trans(1,3); + nums[5]=Trans(2,3); +} + +void FeatureAssociation::runFeatureAssociation(Eigen::Matrix4f Init_state_) { + Matrix4dToRpyxyz(Init_state_,transformCur); + updateTransformation(); +} + +void FeatureAssociation::RpyxyzToMatrix4f(Eigen::Matrix4f &Trans,float nums[]){ + Eigen::Vector3f eulerAngle(nums[0],nums[1],nums[2]); + Eigen::AngleAxisf rollAngle(Eigen::AngleAxisd(eulerAngle(0),Eigen::Vector3d::UnitX())); + Eigen::AngleAxisf pitchAngle(Eigen::AngleAxisd(eulerAngle(1),Eigen::Vector3d::UnitY())); + Eigen::AngleAxisf yawAngle(Eigen::AngleAxisd(eulerAngle(2),Eigen::Vector3d::UnitZ())); + + Eigen::Quaternionf quaternion; + quaternion=rollAngle*pitchAngle*yawAngle; + Eigen::Matrix3f rotate(quaternion); + Eigen::Vector3f tran(nums[3],nums[4],nums[5]); + Trans=Eigen::Matrix4f::Identity(); + Trans.block<3,3>(0,0)=rotate; + Trans.block<3,1>(0,3)=tran; +} + +Eigen::Matrix4f FeatureAssociation::getFinalTransformation() +{ + Eigen::Matrix4f Tran; + RpyxyzToMatrix4f(Tran,transformCur); + return Tran; +} + diff --git a/src/hdl_graph_slam/ford2bag.py b/src/hdl_graph_slam/ford2bag.py new file mode 100755 index 0000000..627bcf8 --- /dev/null +++ b/src/hdl_graph_slam/ford2bag.py @@ -0,0 +1,96 @@ +#!/usr/bin/python +# SPDX-License-Identifier: BSD-2-Clause +import re +import os +import sys +import struct +import numpy +import scipy.io + +import rospy +import rosbag +import progressbar +import sensor_msgs.point_cloud2 as pc2 + +from sensor_msgs.msg import NavSatFix, NavSatStatus, PointCloud2 +from geographic_msgs.msg import GeoPointStamped + + +def gps2navsat(filename, bag): + with open(filename, 'r') as file: + try: + while True: + data = struct.unpack('qddd', file.read(8*4)) + time = data[0] + local = data[1:] + lat_lon_el_theta = struct.unpack('dddd', file.read(8*4)) + gps_cov = struct.unpack('dddddddddddddddd', file.read(8*16)) + + if abs(lat_lon_el_theta[0]) < 1e-1: + continue + + navsat = NavSatFix() + navsat.header.frame_id = 'gps' + navsat.header.stamp = rospy.Time.from_sec(time * 1e-6) + navsat.status.status = NavSatStatus.STATUS_FIX + navsat.status.service = NavSatStatus.SERVICE_GPS + + navsat.latitude = lat_lon_el_theta[0] + navsat.longitude = lat_lon_el_theta[1] + navsat.altitude = lat_lon_el_theta[2] + + navsat.position_covariance = numpy.array(gps_cov).reshape(4, 4)[:3, :3].flatten().tolist() + navsat.position_covariance_type = NavSatFix.COVARIANCE_TYPE_KNOWN + + bag.write('/gps/fix', navsat, navsat.header.stamp) + + geopoint = GeoPointStamped() + geopoint.header = navsat.header + geopoint.position.latitude = lat_lon_el_theta[0] + geopoint.position.longitude = lat_lon_el_theta[1] + geopoint.position.altitude = lat_lon_el_theta[2] + + bag.write('/gps/geopoint', geopoint, geopoint.header.stamp) + + except: + print 'done' + + +def mat2pointcloud(filename): + m = scipy.io.loadmat(filename) + scan = numpy.transpose(m['SCAN']['XYZ'][0][0]).astype(numpy.float32) + stamp = m['SCAN']['timestamp_laser'][0][0][0][0] * 1e-6 + + cloud = PointCloud2() + cloud.header.stamp = rospy.Time.from_sec(stamp) + cloud.header.frame_id = 'velodyne' + cloud = pc2.create_cloud_xyz32(cloud.header, scan) + return cloud + + +def main(): + if len(sys.argv) < 2: + print 'usage: ford2bag.py src_dirname output_filename' + + output_filename = sys.argv[1] + + rospy.init_node('bag') + filenames = sorted(['SCANS/' + x for x in os.listdir('SCANS') if re.match('Scan[0-9]*\.mat', x)]) + print filenames + + progress = progressbar.ProgressBar(max_value=len(filenames)) + pub = rospy.Publisher('/velodyne_points', PointCloud2, queue_size=32) + with rosbag.Bag(output_filename, 'w') as bag: + gps2navsat('GPS.log', bag) + for i, filename in enumerate(filenames): + if rospy.is_shutdown(): + break + progress.update(i) + cloud = mat2pointcloud(filename) + if pub.get_num_connections(): + pub.publish(cloud) + bag.write('/velodyne_points', cloud, cloud.header.stamp) + + +if __name__ == '__main__': + main() diff --git a/src/hdl_graph_slam/graph_slam.cpp b/src/hdl_graph_slam/graph_slam.cpp new file mode 100644 index 0000000..8d94c5b --- /dev/null +++ b/src/hdl_graph_slam/graph_slam.cpp @@ -0,0 +1,351 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +G2O_USE_OPTIMIZATION_LIBRARY(pcg) +G2O_USE_OPTIMIZATION_LIBRARY(cholmod) // be aware of that cholmod brings GPL dependency +G2O_USE_OPTIMIZATION_LIBRARY(csparse) // be aware of that csparse brings LGPL unless it is dynamically linked + +namespace g2o { +G2O_REGISTER_TYPE(EDGE_SE3_PLANE, EdgeSE3Plane) +G2O_REGISTER_TYPE(EDGE_SE3_PRIORXY, EdgeSE3PriorXY) +G2O_REGISTER_TYPE(EDGE_SE3_PRIORXYZ, EdgeSE3PriorXYZ) +G2O_REGISTER_TYPE(EDGE_SE3_PRIORVEC, EdgeSE3PriorVec) +G2O_REGISTER_TYPE(EDGE_SE3_PRIORQUAT, EdgeSE3PriorQuat) +G2O_REGISTER_TYPE(EDGE_PLANE_PRIOR_NORMAL, EdgePlanePriorNormal) +G2O_REGISTER_TYPE(EDGE_PLANE_PRIOR_DISTANCE, EdgePlanePriorDistance) +G2O_REGISTER_TYPE(EDGE_PLANE_IDENTITY, EdgePlaneIdentity) +G2O_REGISTER_TYPE(EDGE_PLANE_PARALLEL, EdgePlaneParallel) +G2O_REGISTER_TYPE(EDGE_PLANE_PAERPENDICULAR, EdgePlanePerpendicular) +} // namespace g2o + +namespace hdl_graph_slam { + +/** + * @brief constructor + */ +GraphSLAM::GraphSLAM(const std::string& solver_type) { + graph.reset(new g2o::SparseOptimizer()); + g2o::SparseOptimizer* graph = dynamic_cast(this->graph.get()); + + std::cout << "construct solver: " << solver_type << std::endl; + g2o::OptimizationAlgorithmFactory* solver_factory = g2o::OptimizationAlgorithmFactory::instance(); + g2o::OptimizationAlgorithmProperty solver_property; + g2o::OptimizationAlgorithm* solver = solver_factory->construct(solver_type, solver_property); + graph->setAlgorithm(solver); + + if(!graph->solver()) { + std::cerr << std::endl; + std::cerr << "error : failed to allocate solver!!" << std::endl; + solver_factory->listSolvers(std::cerr); + std::cerr << "-------------" << std::endl; + std::cin.ignore(1); + return; + } + std::cout << "done" << std::endl; + + robust_kernel_factory = g2o::RobustKernelFactory::instance(); +} + +/** + * @brief destructor + */ +GraphSLAM::~GraphSLAM() { + graph.reset(); +} + +void GraphSLAM::set_solver(const std::string& solver_type) { + g2o::SparseOptimizer* graph = dynamic_cast(this->graph.get()); + + std::cout << "construct solver: " << solver_type << std::endl; + g2o::OptimizationAlgorithmFactory* solver_factory = g2o::OptimizationAlgorithmFactory::instance(); + g2o::OptimizationAlgorithmProperty solver_property; + g2o::OptimizationAlgorithm* solver = solver_factory->construct(solver_type, solver_property); + graph->setAlgorithm(solver); + + if(!graph->solver()) { + std::cerr << std::endl; + std::cerr << "error : failed to allocate solver!!" << std::endl; + solver_factory->listSolvers(std::cerr); + std::cerr << "-------------" << std::endl; + std::cin.ignore(1); + return; + } + std::cout << "done" << std::endl; +} + +int GraphSLAM::num_vertices() const { + return graph->vertices().size(); +} +int GraphSLAM::num_edges() const { + return graph->edges().size(); +} + +g2o::VertexSE3* GraphSLAM::add_se3_node(const Eigen::Isometry3d& pose) { + g2o::VertexSE3* vertex(new g2o::VertexSE3()); + vertex->setId(static_cast(graph->vertices().size())); + vertex->setEstimate(pose); + graph->addVertex(vertex); + + return vertex; +} + +g2o::VertexPlane* GraphSLAM::add_plane_node(const Eigen::Vector4d& plane_coeffs) { + g2o::VertexPlane* vertex(new g2o::VertexPlane()); + vertex->setId(static_cast(graph->vertices().size())); + vertex->setEstimate(plane_coeffs); + graph->addVertex(vertex); + + return vertex; +} + +g2o::VertexPointXYZ* GraphSLAM::add_point_xyz_node(const Eigen::Vector3d& xyz) { + g2o::VertexPointXYZ* vertex(new g2o::VertexPointXYZ()); + vertex->setId(static_cast(graph->vertices().size())); + vertex->setEstimate(xyz); + graph->addVertex(vertex); + + return vertex; +} + +g2o::EdgeSE3* GraphSLAM::add_se3_edge(g2o::VertexSE3* v1, g2o::VertexSE3* v2, const Eigen::Isometry3d& relative_pose, const Eigen::MatrixXd& information_matrix) { + g2o::EdgeSE3* edge(new g2o::EdgeSE3()); + edge->setMeasurement(relative_pose); + edge->setInformation(information_matrix); + edge->vertices()[0] = v1; + edge->vertices()[1] = v2; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgeSE3Plane* GraphSLAM::add_se3_plane_edge(g2o::VertexSE3* v_se3, g2o::VertexPlane* v_plane, const Eigen::Vector4d& plane_coeffs, const Eigen::MatrixXd& information_matrix) { + g2o::EdgeSE3Plane* edge(new g2o::EdgeSE3Plane()); + edge->setMeasurement(plane_coeffs); + edge->setInformation(information_matrix); + edge->vertices()[0] = v_se3; + edge->vertices()[1] = v_plane; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgeSE3PointXYZ* GraphSLAM::add_se3_point_xyz_edge(g2o::VertexSE3* v_se3, g2o::VertexPointXYZ* v_xyz, const Eigen::Vector3d& xyz, const Eigen::MatrixXd& information_matrix) { + g2o::EdgeSE3PointXYZ* edge(new g2o::EdgeSE3PointXYZ()); + edge->setMeasurement(xyz); + edge->setInformation(information_matrix); + edge->vertices()[0] = v_se3; + edge->vertices()[1] = v_xyz; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgePlanePriorNormal* GraphSLAM::add_plane_normal_prior_edge(g2o::VertexPlane* v, const Eigen::Vector3d& normal, const Eigen::MatrixXd& information_matrix) { + g2o::EdgePlanePriorNormal* edge(new g2o::EdgePlanePriorNormal()); + edge->setMeasurement(normal); + edge->setInformation(information_matrix); + edge->vertices()[0] = v; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgePlanePriorDistance* GraphSLAM::add_plane_distance_prior_edge(g2o::VertexPlane* v, double distance, const Eigen::MatrixXd& information_matrix) { + g2o::EdgePlanePriorDistance* edge(new g2o::EdgePlanePriorDistance()); + edge->setMeasurement(distance); + edge->setInformation(information_matrix); + edge->vertices()[0] = v; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgeSE3PriorXY* GraphSLAM::add_se3_prior_xy_edge(g2o::VertexSE3* v_se3, const Eigen::Vector2d& xy, const Eigen::MatrixXd& information_matrix) { + g2o::EdgeSE3PriorXY* edge(new g2o::EdgeSE3PriorXY()); + edge->setMeasurement(xy); + edge->setInformation(information_matrix); + edge->vertices()[0] = v_se3; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgeSE3PriorXYZ* GraphSLAM::add_se3_prior_xyz_edge(g2o::VertexSE3* v_se3, const Eigen::Vector3d& xyz, const Eigen::MatrixXd& information_matrix) { + g2o::EdgeSE3PriorXYZ* edge(new g2o::EdgeSE3PriorXYZ()); + edge->setMeasurement(xyz); + edge->setInformation(information_matrix); + edge->vertices()[0] = v_se3; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgeSE3PriorVec* GraphSLAM::add_se3_prior_vec_edge(g2o::VertexSE3* v_se3, const Eigen::Vector3d& direction, const Eigen::Vector3d& measurement, const Eigen::MatrixXd& information_matrix) { + Eigen::Matrix m; + m.head<3>() = direction; + m.tail<3>() = measurement; + + g2o::EdgeSE3PriorVec* edge(new g2o::EdgeSE3PriorVec()); + edge->setMeasurement(m); + edge->setInformation(information_matrix); + edge->vertices()[0] = v_se3; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgeSE3PriorQuat* GraphSLAM::add_se3_prior_quat_edge(g2o::VertexSE3* v_se3, const Eigen::Quaterniond& quat, const Eigen::MatrixXd& information_matrix) { + g2o::EdgeSE3PriorQuat* edge(new g2o::EdgeSE3PriorQuat()); + edge->setMeasurement(quat); + edge->setInformation(information_matrix); + edge->vertices()[0] = v_se3; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgePlane* GraphSLAM::add_plane_edge(g2o::VertexPlane* v_plane1, g2o::VertexPlane* v_plane2, const Eigen::Vector4d& measurement, const Eigen::Matrix4d& information) { + g2o::EdgePlane* edge(new g2o::EdgePlane()); + edge->setMeasurement(measurement); + edge->setInformation(information); + edge->vertices()[0] = v_plane1; + edge->vertices()[1] = v_plane2; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgePlaneIdentity* GraphSLAM::add_plane_identity_edge(g2o::VertexPlane* v_plane1, g2o::VertexPlane* v_plane2, const Eigen::Vector4d& measurement, const Eigen::Matrix4d& information) { + g2o::EdgePlaneIdentity* edge(new g2o::EdgePlaneIdentity()); + edge->setMeasurement(measurement); + edge->setInformation(information); + edge->vertices()[0] = v_plane1; + edge->vertices()[1] = v_plane2; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgePlaneParallel* GraphSLAM::add_plane_parallel_edge(g2o::VertexPlane* v_plane1, g2o::VertexPlane* v_plane2, const Eigen::Vector3d& measurement, const Eigen::Matrix3d& information) { + g2o::EdgePlaneParallel* edge(new g2o::EdgePlaneParallel()); + edge->setMeasurement(measurement); + edge->setInformation(information); + edge->vertices()[0] = v_plane1; + edge->vertices()[1] = v_plane2; + graph->addEdge(edge); + + return edge; +} + +g2o::EdgePlanePerpendicular* GraphSLAM::add_plane_perpendicular_edge(g2o::VertexPlane* v_plane1, g2o::VertexPlane* v_plane2, const Eigen::Vector3d& measurement, const Eigen::MatrixXd& information) { + g2o::EdgePlanePerpendicular* edge(new g2o::EdgePlanePerpendicular()); + edge->setMeasurement(measurement); + edge->setInformation(information); + edge->vertices()[0] = v_plane1; + edge->vertices()[1] = v_plane2; + graph->addEdge(edge); + + return edge; +} + +void GraphSLAM::add_robust_kernel(g2o::HyperGraph::Edge* edge, const std::string& kernel_type, double kernel_size) { + if(kernel_type == "NONE") { + return; + } + + g2o::RobustKernel* kernel = robust_kernel_factory->construct(kernel_type); + if(kernel == nullptr) { + std::cerr << "warning : invalid robust kernel type: " << kernel_type << std::endl; + return; + } + + kernel->setDelta(kernel_size); + + g2o::OptimizableGraph::Edge* edge_ = dynamic_cast(edge); + edge_->setRobustKernel(kernel); +} + +int GraphSLAM::optimize(int num_iterations) { + g2o::SparseOptimizer* graph = dynamic_cast(this->graph.get()); + if(graph->edges().size() < 10) { + return -1; + } + + //std::cout << std::endl; + //std::cout << "--- pose graph optimization ---" << std::endl; + //std::cout << "nodes: " << graph->vertices().size() << " edges: " << graph->edges().size() << std::endl; + //std::cout << "optimizing... " << std::flush; + + //std::cout << "init" << std::endl; + graph->initializeOptimization(); + graph->setVerbose(false); + + //std::cout << "chi2" << std::endl; + double chi2 = graph->chi2(); + + //std::cout << "optimize!!" << std::endl; + auto t1 = ros::WallTime::now(); + int iterations = graph->optimize(num_iterations); + + auto t2 = ros::WallTime::now(); + //std::cout << "done" << std::endl; + //std::cout << "iterations: " << iterations << " / " << num_iterations << std::endl; + //std::cout << "chi2: (before)" << chi2 << " -> (after)" << graph->chi2() << std::endl; + //std::cout << "time: " << boost::format("%.3f") % (t2 - t1).toSec() << "[sec]" << std::endl; + + return iterations; +} + +void GraphSLAM::save(const std::string& filename) { + g2o::SparseOptimizer* graph = dynamic_cast(this->graph.get()); + + std::ofstream ofs(filename); + graph->save(ofs); + + g2o::save_robust_kernels(filename + ".kernels", graph); +} + +bool GraphSLAM::load(const std::string& filename) { + std::cout << "loading pose graph..." << std::endl; + g2o::SparseOptimizer* graph = dynamic_cast(this->graph.get()); + + std::ifstream ifs(filename); + if(!graph->load(ifs)) { + return false; + } + + std::cout << "nodes : " << graph->vertices().size() << std::endl; + std::cout << "edges : " << graph->edges().size() << std::endl; + + if(!g2o::load_robust_kernels(filename + ".kernels", graph)) { + return false; + } + + return true; +} + +} // namespace hdl_graph_slam diff --git a/src/hdl_graph_slam/imageProjection.cpp b/src/hdl_graph_slam/imageProjection.cpp new file mode 100644 index 0000000..70f3244 --- /dev/null +++ b/src/hdl_graph_slam/imageProjection.cpp @@ -0,0 +1,385 @@ +#include "hdl_graph_slam/imageProjection.h" +#include +#include + +using namespace hdl_graph_slam; +ImageProjection::ImageProjection() { + + nanPoint.x = std::numeric_limits::quiet_NaN(); + nanPoint.y = std::numeric_limits::quiet_NaN(); + nanPoint.z = std::numeric_limits::quiet_NaN(); + nanPoint.intensity = -1; + + allocateMemory(); + resetParameters(); +} + +void ImageProjection::allocateMemory() { + if(Time_<10) + last_distr_.setOnes(); + sigma_smooth_predict_ = 6; + laserCloudIn.reset(new pcl::PointCloud()); + laserCloudIn_last_.reset(new pcl::PointCloud()); + LinearPoints_.reset(new pcl::PointCloud()); + LinearPoints_last_.reset(new pcl::PointCloud()); + fullLinearPoints_.reset(new pcl::PointCloud()); + LinearPoints_Ground.reset(new pcl::PointCloud()); + LinearPoints_Less_Ground.reset(new pcl::PointCloud()); + Less_LinearPoints_.reset(new pcl::PointCloud()); + Less_LinearPoints_D_.reset(new pcl::PointCloud()); + + + fullValidPoints_z_.reset(new pcl::PointCloud()); + fullValidPoints_smoothness_.reset(new pcl::PointCloud()); + + fullCloud.reset(new pcl::PointCloud()); + fullInfoCloud.reset(new pcl::PointCloud()); + + groundCloud.reset(new pcl::PointCloud()); + segmentedCloud.reset(new pcl::PointCloud()); + segmentedCloudPure.reset(new pcl::PointCloud()); + outlierCloud.reset(new pcl::PointCloud()); + + fullCloud->points.resize(N_SCAN * Horizon_SCAN); + fullInfoCloud->points.resize(N_SCAN * Horizon_SCAN); + + segMsg.startRingIndex = new int[N_SCAN]; + segMsg.endRingIndex = new int[N_SCAN]; + + segMsg.segmentedCloudGroundFlag = new bool[N_SCAN * Horizon_SCAN]; + segMsg.segmentedCloudColInd = new int[N_SCAN * Horizon_SCAN]; + segMsg.segmentedCloudRange = new float[N_SCAN * Horizon_SCAN]; + + for (int i = 0; i < N_SCAN * Horizon_SCAN; i++) { + segMsg.segmentedCloudGroundFlag[i] = false; + segMsg.segmentedCloudColInd[i] = 0; + segMsg.segmentedCloudRange[i] = 0; + } + + std::pair neighbor; + neighbor.first = -1; + neighbor.second = 0; + neighborIterator.push_back(neighbor); + neighbor.first = 0; + neighbor.second = 1; + neighborIterator.push_back(neighbor); + neighbor.first = 0; + neighbor.second = -1; + neighborIterator.push_back(neighbor); + neighbor.first = 1; + neighbor.second = 0; + neighborIterator.push_back(neighbor); + + allPushedIndX = new uint16_t[N_SCAN * Horizon_SCAN]; + allPushedIndY = new uint16_t[N_SCAN * Horizon_SCAN]; + + queueIndX = new uint16_t[N_SCAN * Horizon_SCAN]; + queueIndY = new uint16_t[N_SCAN * Horizon_SCAN]; + + +} + +void ImageProjection::resetParameters() { + laserCloudIn->clear(); + //laserCloudIn_last_->clear(); + groundCloud->clear(); + segmentedCloud->clear(); + segmentedCloudPure->clear(); + outlierCloud->clear(); + LinearPoints_->clear(); + Less_LinearPoints_->clear(); + Less_LinearPoints_D_->clear(); + fullLinearPoints_->clear(); + LinearPoints_Ground->clear(); + LinearPoints_Less_Ground->clear(); + + fullValidPoints_z_->clear(); + fullValidPoints_smoothness_->clear(); + + rangeMat = cv::Mat(N_SCAN, Horizon_SCAN, CV_32F, cv::Scalar::all(FLT_MAX)); + groundMat = cv::Mat(N_SCAN, Horizon_SCAN, CV_8S, cv::Scalar::all(0)); + labelMat = cv::Mat(N_SCAN, Horizon_SCAN, CV_32S, cv::Scalar::all(0)); + labelCount = 1; + + std::fill(fullCloud->points.begin(), fullCloud->points.end(), nanPoint); + std::fill(fullInfoCloud->points.begin(), fullInfoCloud->points.end(), nanPoint); +} + +void ImageProjection::showTwoPC(pcl::PointCloud::Ptr curr_cloud_gicp, pcl::PointCloud::Ptr ref_cloud) { + + bool visualize = true; + + if (visualize) { + int v1(0); + using pcl::visualization::PCLVisualizer; + using pcl::visualization::PointCloudColorHandlerGenericField; + boost::shared_ptr viewer(new pcl::visualization::PCLVisualizer("3D Viewer")); + viewer->setBackgroundColor(0, 0, 0); + pcl::PointCloud::Ptr curr_cloud_ptr(curr_cloud_gicp);//icp_output + pcl::PointCloud::Ptr ref_cloud_ptr(ref_cloud); + PointCloudColorHandlerGenericField fildColor(curr_cloud_ptr, "z"); + pcl::visualization::PointCloudColorHandlerCustom single_color0(curr_cloud_ptr, 0, 255, 0); // green + + viewer->addPointCloud(curr_cloud_ptr, single_color0, "sample cloud"); + pcl::visualization::PointCloudColorHandlerCustom single_color1(ref_cloud_ptr, 255, 0, 0); // red + viewer->addPointCloud(ref_cloud_ptr, single_color1, "sample cloud1"); + viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 1, "sample cloud"); + //viewer->createViewPort(0,0,1,1,v1); + viewer->addCoordinateSystem(1.0); + viewer->initCameraParameters(); + + while (!viewer->wasStopped()) { + viewer->spinOnce(100); + //boost::this_thread::sleep(boost::posix_time::microseconds(100000)); + } + } +} + + +void ImageProjection::cloudHandler(pcl::PointCloud::Ptr input_cloud) { + + clock_t start0 = clock(); + for (int i = 0; i < input_cloud->size(); i++) + laserCloudIn->push_back(input_cloud->points[i]); + + clock_t start = clock(); + //std::cout<<"get pc"<<(double)(start-start0)/CLOCKS_PER_SEC<points[0].y, laserCloudIn->points[0].x); + segMsg.endOrientation = -atan2(laserCloudIn->points[laserCloudIn->points.size() - 1].y, + laserCloudIn->points[laserCloudIn->points.size() - 1].x) + 2 * M_PI; + if (segMsg.endOrientation - segMsg.startOrientation > 3 * M_PI) { + segMsg.endOrientation -= 2 * M_PI; + } else if (segMsg.endOrientation - segMsg.startOrientation < M_PI) + segMsg.endOrientation += 2 * M_PI; + segMsg.orientationDiff = segMsg.endOrientation - segMsg.startOrientation; +} + +void ImageProjection::projectPointCloud() { + // range image projection + float verticalAngle, horizonAngle, range; + size_t rowIdn, columnIdn, index, cloudSize; + PointType thisPoint; + + cloudSize = laserCloudIn->points.size(); + + for (size_t i = 0; i < cloudSize; ++i) { + + thisPoint.x = laserCloudIn->points[i].x; + thisPoint.y = laserCloudIn->points[i].y; + thisPoint.z = laserCloudIn->points[i].z; + // find the row and column index in the image for this point + verticalAngle = atan2(thisPoint.z, sqrt(thisPoint.x * thisPoint.x + thisPoint.y * thisPoint.y)) * 180 / M_PI; + rowIdn = (verticalAngle + ang_bottom) / ang_res_y; + if (rowIdn < 0 || rowIdn >= N_SCAN) + continue; + + horizonAngle = atan2(thisPoint.x, thisPoint.y) * 180 / M_PI; + + columnIdn = -round((horizonAngle - 90.0) / ang_res_x) + Horizon_SCAN / 2; + if (columnIdn >= Horizon_SCAN) + columnIdn -= Horizon_SCAN; + + if (columnIdn < 0 || columnIdn >= Horizon_SCAN) + continue; + + range = sqrt(thisPoint.x * thisPoint.x + thisPoint.y * thisPoint.y + thisPoint.z * thisPoint.z); + if (range < 0.5) + continue; + + rangeMat.at(rowIdn, columnIdn) = range; + + thisPoint.intensity = (float) rowIdn + (float) columnIdn / 10000.0; + + index = columnIdn + rowIdn * Horizon_SCAN; + fullCloud->points[index] = thisPoint; + fullInfoCloud->points[index] = thisPoint; + fullInfoCloud->points[index].intensity = range; // the corresponding range of a point is saved as "intensity" + } + +} + + +void ImageProjection::extractLinearPoints() { + + use_prob_smooth_ = true; + size_t upperInd, currInd, leftInd_1, rightInd_1, leftInd_2, rightInd_2, leftInd_3, + rightInd_3, lowInd,rightupperInd_1,leftupperInd_1; + int point_no = 0, point_no_g(0), point_no_less_g(0), point_no_less = 0,point_no_full(0); + + plane_no_ = 0; + + + for (size_t j = 3; j < Horizon_SCAN - 4; j = j+1) { + for (size_t i = 1; i < N_SCAN - 1; ++i) { + lowInd = j + (i - 1) * Horizon_SCAN; + currInd = j + i * Horizon_SCAN; + upperInd = j + (i + 1) * Horizon_SCAN; + leftInd_1 = j - 1 + i * Horizon_SCAN; + leftupperInd_1 = j - 1 + (i+1) * Horizon_SCAN; + rightInd_1 = j + 1 + i * Horizon_SCAN; + rightupperInd_1 = j + 1 + (i+1) * Horizon_SCAN; + leftInd_2 = j - 2 + i * Horizon_SCAN; + rightInd_2 = j + 2 + i * Horizon_SCAN; + leftInd_3 = j - 3 + i * Horizon_SCAN; + rightInd_3 = j + 3 + i * Horizon_SCAN; + + if (fullInfoCloud->points[currInd].intensity == -1 || + fullInfoCloud->points[upperInd].intensity == -1 || + fullInfoCloud->points[leftInd_1].intensity == -1 || + fullInfoCloud->points[rightInd_1].intensity == -1 + ) { + continue; + } + + double smoothness = + fullInfoCloud->points[leftInd_1].intensity + fullInfoCloud->points[leftInd_2].intensity + + fullInfoCloud->points[leftInd_3].intensity - 6 * fullInfoCloud->points[currInd].intensity + + fullInfoCloud->points[rightInd_1].intensity + fullInfoCloud->points[rightInd_2].intensity + + fullInfoCloud->points[rightInd_3].intensity; + PointType upper_normals0; + upper_normals0.x = fullInfoCloud->points[upperInd].x - fullInfoCloud->points[currInd].x; + upper_normals0.y = fullInfoCloud->points[upperInd].y - fullInfoCloud->points[currInd].y; + upper_normals0.z = fullInfoCloud->points[upperInd].z - fullInfoCloud->points[currInd].z; + upper_normals0.intensity = 0; + + float norm = sqrt( + upper_normals0.x * upper_normals0.x + upper_normals0.y * upper_normals0.y + upper_normals0.z * + upper_normals0.z); + upper_normals0.z = upper_normals0.z / norm; + smoothness = abs(smoothness) / 6; + + + float z_thre = 0.9; + if (abs(upper_normals0.z) > z_thre and groundMat.at(i, j) > -1) // + { + PointType p = fullInfoCloud->points[currInd]; + p.intensity = smoothness; + Less_LinearPoints_->points.push_back(p); + point_no_less++; + p.intensity = smoothness; + if(smoothness > 1 and currInd%2 == 0) //and currInd%2 == 0 + { + LinearPoints_->points.push_back(p); + ++point_no; + } + } + + if(abs(upper_normals0.z)<0.3 ) // and groundMat.at(i, j) == 1 + { + if (currInd % 2 == 0) { + LinearPoints_Less_Ground->points.push_back(fullInfoCloud->points[currInd]); + ++point_no_less_g; + } + if (currInd % 4 == 0) { + LinearPoints_Ground->points.push_back(fullInfoCloud->points[currInd]); + ++point_no_g; + } + } + } + } +/* + for (size_t j = 0; j < Horizon_SCAN; ++j) { + for (size_t i = 0; i < groundScanInd; ++i) { + size_t currInd = j + (i) * Horizon_SCAN; + if(groundMat.at(i, j) == 1) { + if (currInd % 4 == 0) { + LinearPoints_Less_Ground->points.push_back(fullInfoCloud->points[currInd]); + ++point_no_less_g; + } + if (currInd % 16 == 0) { + LinearPoints_Ground->points.push_back(fullInfoCloud->points[currInd]); + ++point_no_g; + } + } + } + }*/ + + //smooth_f.close(); + LinearPoints_->width = point_no; + LinearPoints_->height = 1; + LinearPoints_->is_dense = false; + LinearPoints_->points.resize(LinearPoints_->width * LinearPoints_->height); + // cout<<"LinearPoints_->size()"<<" "<size()<<" "<size()<width = point_no_less; + Less_LinearPoints_->height = 1; + Less_LinearPoints_->is_dense = false; + Less_LinearPoints_->points.resize(Less_LinearPoints_->width * Less_LinearPoints_->height); + //cout<<"Less_LinearPoints_->size()"<<" "<size()<width = point_no_g; + LinearPoints_Ground->height = 1; + LinearPoints_Ground->is_dense = false; + LinearPoints_Ground->points.resize(LinearPoints_Ground->width * LinearPoints_Ground->height); + // cout<<"LinearPoints_Ground->size()"<<" "<size()<width = point_no_less_g; + LinearPoints_Less_Ground->height = 1; + LinearPoints_Less_Ground->is_dense = false; + LinearPoints_Less_Ground->points.resize(LinearPoints_Less_Ground->width * LinearPoints_Less_Ground->height); + //cout<<"LinearPoints_Less_Ground->size()"<<" "<size()<points[lowerInd].intensity == -1 || + fullCloud->points[upperInd].intensity == -1) { + // no info to check, invalid points + groundMat.at(i, j) = -1; + continue; + } + + diffX = fullCloud->points[upperInd].x - fullCloud->points[lowerInd].x; + diffY = fullCloud->points[upperInd].y - fullCloud->points[lowerInd].y; + diffZ = fullCloud->points[upperInd].z - fullCloud->points[lowerInd].z; + + angle = atan2(diffZ, sqrt(diffX * diffX + diffY * diffY)) * 180 / M_PI; + + if (abs(angle - sensorMountAngle) <= 10) { + groundMat.at(i, j) = 1; + groundMat.at(i + 1, j) = 1; + + } + } + } + // extract ground cloud (groundMat == 1) + // mark entry that doesn't need to label (ground and invalid point) for segmentation + // note that ground remove is from 0~N_SCAN-1, need rangeMat for mark label matrix for the 16th scan + for (size_t i = 0; i < N_SCAN; ++i) { + for (size_t j = 0; j < Horizon_SCAN; ++j) { + if (groundMat.at(i, j) == 1 || rangeMat.at(i, j) == FLT_MAX) { + labelMat.at(i, j) = -1; + } + } + } +} diff --git a/src/hdl_graph_slam/information_matrix_calculator.cpp b/src/hdl_graph_slam/information_matrix_calculator.cpp new file mode 100644 index 0000000..f3a7e1a --- /dev/null +++ b/src/hdl_graph_slam/information_matrix_calculator.cpp @@ -0,0 +1,82 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include + +#include +#include + +namespace hdl_graph_slam { + +InformationMatrixCalculator::InformationMatrixCalculator(ros::NodeHandle& nh) { + use_const_inf_matrix = nh.param("use_const_inf_matrix", false); + const_stddev_x = nh.param("const_stddev_x", 0.5); + const_stddev_q = nh.param("const_stddev_q", 0.1); + + var_gain_a = nh.param("var_gain_a", 20.0); + min_stddev_x = nh.param("min_stddev_x", 0.1); + max_stddev_x = nh.param("max_stddev_x", 5.0); + min_stddev_q = nh.param("min_stddev_q", 0.05); + max_stddev_q = nh.param("max_stddev_q", 0.2); + fitness_score_thresh = nh.param("fitness_score_thresh", 0.5); +} + +InformationMatrixCalculator::~InformationMatrixCalculator() {} + +Eigen::MatrixXd InformationMatrixCalculator::calc_information_matrix(const pcl::PointCloud::ConstPtr& cloud1, const pcl::PointCloud::ConstPtr& cloud2, const Eigen::Isometry3d& relpose) const { + if(use_const_inf_matrix) { + Eigen::MatrixXd inf = Eigen::MatrixXd::Identity(6, 6); + inf.topLeftCorner(3, 3).array() /= const_stddev_x; + inf.bottomRightCorner(3, 3).array() /= const_stddev_q; + return inf; + } + + double fitness_score = calc_fitness_score(cloud1, cloud2, relpose); + + double min_var_x = std::pow(min_stddev_x, 2); + double max_var_x = std::pow(max_stddev_x, 2); + double min_var_q = std::pow(min_stddev_q, 2); + double max_var_q = std::pow(max_stddev_q, 2); + + float w_x = weight(var_gain_a, fitness_score_thresh, min_var_x, max_var_x, fitness_score); + float w_q = weight(var_gain_a, fitness_score_thresh, min_var_q, max_var_q, fitness_score); + + Eigen::MatrixXd inf = Eigen::MatrixXd::Identity(6, 6); + inf.topLeftCorner(3, 3).array() /= w_x; + inf.bottomRightCorner(3, 3).array() /= w_q; + return inf; +} + +double InformationMatrixCalculator::calc_fitness_score(const pcl::PointCloud::ConstPtr& cloud1, const pcl::PointCloud::ConstPtr& cloud2, const Eigen::Isometry3d& relpose, double max_range) { + pcl::search::KdTree::Ptr tree_(new pcl::search::KdTree()); + tree_->setInputCloud(cloud1); + + double fitness_score = 0.0; + + // Transform the input dataset using the final transformation + pcl::PointCloud input_transformed; + pcl::transformPointCloud(*cloud2, input_transformed, relpose.cast()); + + std::vector nn_indices(1); + std::vector nn_dists(1); + + // For each point in the source dataset + int nr = 0; + for(size_t i = 0; i < input_transformed.points.size(); ++i) { + // Find its nearest neighbor in the target + tree_->nearestKSearch(input_transformed.points[i], 1, nn_indices, nn_dists); + + // Deal with occlusions (incomplete targets) + if(nn_dists[0] <= max_range) { + // Add to the fitness score + fitness_score += nn_dists[0]; + nr++; + } + } + + if(nr > 0) + return (fitness_score / nr); + else + return (std::numeric_limits::max()); +} + +} // namespace hdl_graph_slam diff --git a/src/hdl_graph_slam/keyframe.cpp b/src/hdl_graph_slam/keyframe.cpp new file mode 100644 index 0000000..99edbae --- /dev/null +++ b/src/hdl_graph_slam/keyframe.cpp @@ -0,0 +1,161 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include + +#include + +#include +#include +#include + +namespace hdl_graph_slam { + +KeyFrame::KeyFrame(const ros::Time& stamp, const Eigen::Isometry3d& odom, double accum_distance, const pcl::PointCloud::ConstPtr& cloud) : stamp(stamp), odom(odom), accum_distance(accum_distance), cloud(cloud), node(nullptr) {} + +KeyFrame::KeyFrame(const std::string& directory, g2o::HyperGraph* graph) : stamp(), odom(Eigen::Isometry3d::Identity()), accum_distance(-1), cloud(nullptr), node(nullptr) { + load(directory, graph); +} + +KeyFrame::~KeyFrame() {} + +void KeyFrame::save(const std::string& directory) { + if(!boost::filesystem::is_directory(directory)) { + boost::filesystem::create_directory(directory); + } + + std::ofstream ofs(directory + "/data"); + ofs << "stamp " << stamp.sec << " " << stamp.nsec << "\n"; + + ofs << "estimate\n"; + ofs << node->estimate().matrix() << "\n"; + + ofs << "odom\n"; + ofs << odom.matrix() << "\n"; + + ofs << "accum_distance " << accum_distance << "\n"; + + if(floor_coeffs) { + ofs << "floor_coeffs " << floor_coeffs->transpose() << "\n"; + } + + if(utm_coord) { + ofs << "utm_coord " << utm_coord->transpose() << "\n"; + } + + if(acceleration) { + ofs << "acceleration " << acceleration->transpose() << "\n"; + } + + if(orientation) { + ofs << "orientation " << orientation->w() << " " << orientation->x() << " " << orientation->y() << " " << orientation->z() << "\n"; + } + + if(node) { + ofs << "id " << node->id() << "\n"; + } + + pcl::io::savePCDFileBinary(directory + "/cloud.pcd", *cloud); +} + +bool KeyFrame::load(const std::string& directory, g2o::HyperGraph* graph) { + std::ifstream ifs(directory + "/data"); + if(!ifs) { + return false; + } + + long node_id = -1; + boost::optional estimate; + + while(!ifs.eof()) { + std::string token; + ifs >> token; + + if(token == "stamp") { + ifs >> stamp.sec >> stamp.nsec; + } else if(token == "estimate") { + Eigen::Matrix4d mat; + for(int i = 0; i < 4; i++) { + for(int j = 0; j < 4; j++) { + ifs >> mat(i, j); + } + } + estimate = Eigen::Isometry3d::Identity(); + estimate->linear() = mat.block<3, 3>(0, 0); + estimate->translation() = mat.block<3, 1>(0, 3); + } else if(token == "odom") { + Eigen::Matrix4d odom_mat = Eigen::Matrix4d::Identity(); + for(int i = 0; i < 4; i++) { + for(int j = 0; j < 4; j++) { + ifs >> odom_mat(i, j); + } + } + + odom.setIdentity(); + odom.linear() = odom_mat.block<3, 3>(0, 0); + odom.translation() = odom_mat.block<3, 1>(0, 3); + } else if(token == "accum_distance") { + ifs >> accum_distance; + } else if(token == "floor_coeffs") { + Eigen::Vector4d coeffs; + ifs >> coeffs[0] >> coeffs[1] >> coeffs[2] >> coeffs[3]; + floor_coeffs = coeffs; + } else if(token == "utm_coord") { + Eigen::Vector3d coord; + ifs >> coord[0] >> coord[1] >> coord[2]; + utm_coord = coord; + } else if(token == "acceleration") { + Eigen::Vector3d acc; + ifs >> acc[0] >> acc[1] >> acc[2]; + acceleration = acc; + } else if(token == "orientation") { + Eigen::Quaterniond quat; + ifs >> quat.w() >> quat.x() >> quat.y() >> quat.z(); + orientation = quat; + } else if(token == "id") { + ifs >> node_id; + } + } + + if(node_id < 0) { + ROS_ERROR_STREAM("invalid node id!!"); + ROS_ERROR_STREAM(directory); + return false; + } + + if(graph->vertices().find(node_id) == graph->vertices().end()) { + ROS_ERROR_STREAM("vertex ID=" << node_id << " does not exist!!"); + return false; + } + + node = dynamic_cast(graph->vertices()[node_id]); + if(node == nullptr) { + ROS_ERROR_STREAM("failed to downcast!!"); + return false; + } + + if(estimate) { + node->setEstimate(*estimate); + } + + pcl::PointCloud::Ptr cloud_(new pcl::PointCloud()); + pcl::io::loadPCDFile(directory + "/cloud.pcd", *cloud_); + cloud = cloud_; + + return true; +} + +long KeyFrame::id() const { + return node->id(); +} + +Eigen::Isometry3d KeyFrame::estimate() const { + return node->estimate(); +} + +KeyFrameSnapshot::KeyFrameSnapshot(const Eigen::Isometry3d& pose, const pcl::PointCloud::ConstPtr& cloud) : pose(pose), cloud(cloud) {} + +KeyFrameSnapshot::KeyFrameSnapshot(const KeyFrame::Ptr& key) : pose(key->node->estimate()), cloud(key->cloud) {} + +KeyFrameSnapshot::~KeyFrameSnapshot() {} + +} // namespace hdl_graph_slam diff --git a/src/hdl_graph_slam/map2odom_publisher.py b/src/hdl_graph_slam/map2odom_publisher.py new file mode 100755 index 0000000..18dcea4 --- /dev/null +++ b/src/hdl_graph_slam/map2odom_publisher.py @@ -0,0 +1,41 @@ +#!/usr/bin/python +# SPDX-License-Identifier: BSD-2-Clause +import tf +import rospy +from geometry_msgs.msg import * + + +class Map2OdomPublisher: + def __init__(self): + self.broadcaster = tf.TransformBroadcaster() + self.subscriber = rospy.Subscriber('/hdl_graph_slam/odom2pub', TransformStamped, self.callback) + + def callback(self, odom_msg): + self.odom_msg = odom_msg + + def spin(self): + if not hasattr(self, 'odom_msg'): + self.broadcaster.sendTransform((0, 0, 0), (0, 0, 0, 1), rospy.Time.now(), 'odom', 'map') + return + + pose = self.odom_msg.transform + pos = (pose.translation.x, pose.translation.y, pose.translation.z) + quat = (pose.rotation.x, pose.rotation.y, pose.rotation.z, pose.rotation.w) + + map_frame_id = self.odom_msg.header.frame_id + odom_frame_id = self.odom_msg.child_frame_id + + self.broadcaster.sendTransform(pos, quat, rospy.Time.now(), odom_frame_id, map_frame_id) + + +def main(): + rospy.init_node('map2odom_publisher') + node = Map2OdomPublisher() + + rate = rospy.Rate(10.0) + while not rospy.is_shutdown(): + node.spin() + rate.sleep() + +if __name__ == '__main__': + main() diff --git a/src/hdl_graph_slam/map_cloud_generator.cpp b/src/hdl_graph_slam/map_cloud_generator.cpp new file mode 100644 index 0000000..cc088de --- /dev/null +++ b/src/hdl_graph_slam/map_cloud_generator.cpp @@ -0,0 +1,55 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include + +#include + +namespace hdl_graph_slam { + +MapCloudGenerator::MapCloudGenerator() {} + +MapCloudGenerator::~MapCloudGenerator() {} + +pcl::PointCloud::Ptr MapCloudGenerator::generate(const std::vector& keyframes, double resolution) const { + if(keyframes.empty()) { + std::cerr << "warning: keyframes empty!!" << std::endl; + return nullptr; + } + + pcl::PointCloud::Ptr cloud(new pcl::PointCloud()); + cloud->reserve(keyframes.front()->cloud->size() * keyframes.size()); + + for(const auto& keyframe : keyframes) { + Eigen::Matrix4f pose = keyframe->pose.matrix().cast(); + //std::cout<<"keyframe cloud for map: "<cloud->points.size()<cloud->points) { + PointT dst_pt; + dst_pt.getVector4fMap() = pose * src_pt.getVector4fMap(); + dst_pt.intensity = src_pt.intensity; + cloud->push_back(dst_pt); + } + } + + cloud->width = cloud->size(); + cloud->height = 1; + cloud->is_dense = false; + + if (resolution <=0.0) + return cloud; // To get unfiltered point cloud with intensity + + pcl::octree::OctreePointCloud octree(resolution); + octree.setInputCloud(cloud); + octree.addPointsFromInputCloud(); + + pcl::PointCloud::Ptr filtered(new pcl::PointCloud()); + octree.getOccupiedVoxelCenters(filtered->points); + + filtered->width = filtered->size(); + filtered->height = 1; + filtered->is_dense = false; + + return filtered; +} + +} // namespace hdl_graph_slam diff --git a/src/hdl_graph_slam/registrations.cpp b/src/hdl_graph_slam/registrations.cpp new file mode 100644 index 0000000..66ebf26 --- /dev/null +++ b/src/hdl_graph_slam/registrations.cpp @@ -0,0 +1,129 @@ +// SPDX-License-Identifier: BSD-2-Clause + +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#ifdef USE_VGICP_CUDA +#include +#endif + +namespace hdl_graph_slam { + +pcl::Registration::Ptr select_registration_method(ros::NodeHandle& pnh) { + using PointT = pcl::PointXYZI; + + // select a registration method (ICP, GICP, NDT) + std::string registration_method = pnh.param("registration_method", "NDT_OMP"); + if(registration_method == "FAST_GICP") { + std::cout << "registration: FAST_GICP" << std::endl; + fast_gicp::FastGICP::Ptr gicp(new fast_gicp::FastGICP()); + gicp->setNumThreads(pnh.param("reg_num_threads", 0)); + gicp->setTransformationEpsilon(pnh.param("reg_transformation_epsilon", 0.01)); + gicp->setMaximumIterations(pnh.param("reg_maximum_iterations", 64)); + gicp->setMaxCorrespondenceDistance(pnh.param("reg_max_correspondence_distance", 2.5)); + gicp->setCorrespondenceRandomness(pnh.param("reg_correspondence_randomness", 20)); + bool use_test = pnh.param("use_prev_test_", false); + std::cout<<"registration->setTest: "<setTest(use_test); + return gicp; + } +#ifdef USE_VGICP_CUDA + else if(registration_method == "FAST_VGICP_CUDA") { + std::cout << "registration: FAST_VGICP_CUDA" << std::endl; + fast_gicp::FastVGICPCuda::Ptr vgicp(new fast_gicp::FastVGICPCuda()); + vgicp->setResolution(pnh.param("reg_resolution", 1.0)); + vgicp->setTransformationEpsilon(pnh.param("reg_transformation_epsilon", 0.01)); + vgicp->setMaximumIterations(pnh.param("reg_maximum_iterations", 64)); + vgicp->setCorrespondenceRandomness(pnh.param("reg_correspondence_randomness", 20)); + return vgicp; + } +#endif + else if(registration_method == "FAST_VGICP") { + std::cout << "registration: FAST_VGICP" << std::endl; + fast_gicp::FastVGICP::Ptr vgicp(new fast_gicp::FastVGICP()); + vgicp->setNumThreads(pnh.param("reg_num_threads", 0)); + vgicp->setResolution(pnh.param("reg_resolution", 1.0)); + vgicp->setTransformationEpsilon(pnh.param("reg_transformation_epsilon", 0.01)); + vgicp->setMaximumIterations(pnh.param("reg_maximum_iterations", 64)); + vgicp->setCorrespondenceRandomness(pnh.param("reg_correspondence_randomness", 20)); + return vgicp; + } else if(registration_method == "ICP") { + std::cout << "registration: ICP" << std::endl; + pcl::IterativeClosestPoint::Ptr icp(new pcl::IterativeClosestPoint()); + icp->setTransformationEpsilon(pnh.param("reg_transformation_epsilon", 0.01)); + icp->setMaximumIterations(pnh.param("reg_maximum_iterations", 64)); + icp->setMaxCorrespondenceDistance(pnh.param("reg_max_correspondence_distance", 2.5)); + icp->setUseReciprocalCorrespondences(pnh.param("reg_use_reciprocal_correspondences", false)); + return icp; + } else if(registration_method.find("GICP") != std::string::npos) { + if(registration_method.find("OMP") == std::string::npos) { + std::cout << "registration: GICP" << std::endl; + pcl::GeneralizedIterativeClosestPoint::Ptr gicp(new pcl::GeneralizedIterativeClosestPoint()); + gicp->setTransformationEpsilon(pnh.param("reg_transformation_epsilon", 0.01)); + gicp->setMaximumIterations(pnh.param("reg_maximum_iterations", 64)); + gicp->setUseReciprocalCorrespondences(pnh.param("reg_use_reciprocal_correspondences", false)); + gicp->setMaxCorrespondenceDistance(pnh.param("reg_max_correspondence_distance", 2.5)); + gicp->setCorrespondenceRandomness(pnh.param("reg_correspondence_randomness", 20)); + gicp->setMaximumOptimizerIterations(pnh.param("reg_max_optimizer_iterations", 20)); + return gicp; + } else { + std::cout << "registration: GICP_OMP" << std::endl; + pclomp::GeneralizedIterativeClosestPoint::Ptr gicp(new pclomp::GeneralizedIterativeClosestPoint()); + gicp->setTransformationEpsilon(pnh.param("reg_transformation_epsilon", 0.01)); + gicp->setMaximumIterations(pnh.param("reg_maximum_iterations", 64)); + gicp->setUseReciprocalCorrespondences(pnh.param("reg_use_reciprocal_correspondences", false)); + gicp->setMaxCorrespondenceDistance(pnh.param("reg_max_correspondence_distance", 2.5)); + gicp->setCorrespondenceRandomness(pnh.param("reg_correspondence_randomness", 20)); + gicp->setMaximumOptimizerIterations(pnh.param("reg_max_optimizer_iterations", 20)); + return gicp; + } + } else { + if(registration_method.find("NDT") == std::string::npos) { + std::cerr << "warning: unknown registration type(" << registration_method << ")" << std::endl; + std::cerr << " : use NDT" << std::endl; + } + + double ndt_resolution = pnh.param("reg_resolution", 0.5); + if(registration_method.find("OMP") == std::string::npos) { + std::cout << "registration: NDT " << ndt_resolution << std::endl; + pcl::NormalDistributionsTransform::Ptr ndt(new pcl::NormalDistributionsTransform()); + ndt->setTransformationEpsilon(pnh.param("reg_transformation_epsilon", 0.01)); + ndt->setMaximumIterations(pnh.param("reg_maximum_iterations", 64)); + ndt->setResolution(ndt_resolution); + return ndt; + } else { + int num_threads = pnh.param("reg_num_threads", 0); + std::string nn_search_method = pnh.param("reg_nn_search_method", "DIRECT7"); + std::cout << "registration: NDT_OMP " << nn_search_method << " " << ndt_resolution << " (" << num_threads << " threads)" << std::endl; + pclomp::NormalDistributionsTransform::Ptr ndt(new pclomp::NormalDistributionsTransform()); + if(num_threads > 0) { + ndt->setNumThreads(num_threads); + } + ndt->setTransformationEpsilon(pnh.param("reg_transformation_epsilon", 0.01)); + ndt->setMaximumIterations(pnh.param("reg_maximum_iterations", 64)); + ndt->setResolution(ndt_resolution); + if(nn_search_method == "KDTREE") { + ndt->setNeighborhoodSearchMethod(pclomp::KDTREE); + } else if(nn_search_method == "DIRECT1") { + ndt->setNeighborhoodSearchMethod(pclomp::DIRECT1); + } else { + ndt->setNeighborhoodSearchMethod(pclomp::DIRECT7); + } + return ndt; + } + } + + return nullptr; +} + +} // namespace hdl_graph_slam diff --git a/srv/DumpGraph.srv b/srv/DumpGraph.srv new file mode 100644 index 0000000..ba7eac2 --- /dev/null +++ b/srv/DumpGraph.srv @@ -0,0 +1,4 @@ + +string destination +--- +bool success diff --git a/srv/SaveMap.srv b/srv/SaveMap.srv new file mode 100644 index 0000000..4fd8b97 --- /dev/null +++ b/srv/SaveMap.srv @@ -0,0 +1,5 @@ +bool utm +float32 resolution +string destination +--- +bool success