diff --git a/README.en.md b/README.en.md
index 397bcfe..3f934f0 100644
--- a/README.en.md
+++ b/README.en.md
@@ -22,53 +22,57 @@ Xiao, K., Ma, L., Tan, S., Cong, Y., Wang, X.: Implementation of UAV Coordinatio
Developers can quickly verify algorithms with XTDrone, such as:
-1. Stereo SLAM
+Stereo SLAM
-2. VIO
+VIO
-3. Dense Reconstruction
+Dense Reconstruction
-4. 2D Laser SLAM
+2D Laser SLAM
-5. 3D Laser SLAM
+3D Laser SLAM
-6. Motion Planning
+2D Motion Planning
-
+
-7. Object Detection and Tracking
+3D Motion Planning
+
+
+
+Object Detection and Tracking
-8. Formation
+Formation
-9. Fixed wing
+Fixed wing
-10. VTOL
+VTOL
-11. UGV
+UGV
-12. USV
+USV
diff --git a/README.md b/README.md
index e422bb1..3664bf9 100644
--- a/README.md
+++ b/README.md
@@ -21,54 +21,58 @@ Xiao, K., Ma, L., Tan, S., Cong, Y., Wang, X.: Implementation of UAV Coordinatio
在这个平台上,开发者可以快速验证算法。如:
-1. 双目SLAM
+双目SLAM
-2. 视觉惯性导航
+视觉惯性导航
-3. 三维稠密重建
+三维稠密重建
-4. 2D激光SLAM
+2D激光SLAM
-5. 3D激光SLAM
+3D激光SLAM
-6. 运动规划
+2D运动规划
-
+
-7. 目标检测与追踪
+3D运动规划
+
+
+
+目标检测与追踪
-8. 多机协同
+多机协同
-9. 固定翼
+固定翼
-10. 复合翼
+复合翼
-11. 无人车
+无人车
-12. 无人船
+无人船
diff --git a/communication/multirotor_communication.py b/communication/multirotor_communication.py
index 0c2a8d3..6843d0b 100644
--- a/communication/multirotor_communication.py
+++ b/communication/multirotor_communication.py
@@ -119,11 +119,15 @@ class Communication:
def cmd_pose_flu_callback(self, msg):
self.coordinate_frame = 9
- self.target_motion = self.construct_target(x=msg.position.x,y=msg.position.y,z=msg.position.z)
+ self.motion_type = 0
+ yaw = self.q2yaw(msg.orientation)
+ self.target_motion = self.construct_target(x=msg.position.x,y=msg.position.y,z=msg.position.z,yaw=yaw)
def cmd_pose_enu_callback(self, msg):
self.coordinate_frame = 1
- self.target_motion = self.construct_target(x=msg.position.x,y=msg.position.y,z=msg.position.z)
+ self.motion_type = 0
+ yaw = self.q2yaw(msg.orientation)
+ self.target_motion = self.construct_target(x=msg.position.x,y=msg.position.y,z=msg.position.z,yaw=yaw)
def cmd_vel_flu_callback(self, msg):
self.hover_state_transition(msg.linear.x, msg.linear.y, msg.linear.z, msg.angular.z)
diff --git a/images/motion_planning.gif b/images/2d_motion_planning.gif
similarity index 100%
rename from images/motion_planning.gif
rename to images/2d_motion_planning.gif
diff --git a/images/3d_motion_planning.gif b/images/3d_motion_planning.gif
new file mode 100755
index 0000000..1e03f71
Binary files /dev/null and b/images/3d_motion_planning.gif differ
diff --git a/motion_planning/launch/2d_motion_planning.launch b/motion_planning/2d/launch/2d_motion_planning.launch
similarity index 100%
rename from motion_planning/launch/2d_motion_planning.launch
rename to motion_planning/2d/launch/2d_motion_planning.launch
diff --git a/motion_planning/launch/move_base.launch b/motion_planning/2d/launch/move_base.launch
similarity index 100%
rename from motion_planning/launch/move_base.launch
rename to motion_planning/2d/launch/move_base.launch
diff --git a/motion_planning/map/indoor3.pgm b/motion_planning/2d/map/indoor3.pgm
similarity index 100%
rename from motion_planning/map/indoor3.pgm
rename to motion_planning/2d/map/indoor3.pgm
diff --git a/motion_planning/map/indoor3.yaml b/motion_planning/2d/map/indoor3.yaml
similarity index 100%
rename from motion_planning/map/indoor3.yaml
rename to motion_planning/2d/map/indoor3.yaml
diff --git a/motion_planning/map/indoor5.pgm b/motion_planning/2d/map/indoor5.pgm
similarity index 100%
rename from motion_planning/map/indoor5.pgm
rename to motion_planning/2d/map/indoor5.pgm
diff --git a/motion_planning/map/indoor5.yaml b/motion_planning/2d/map/indoor5.yaml
similarity index 100%
rename from motion_planning/map/indoor5.yaml
rename to motion_planning/2d/map/indoor5.yaml
diff --git a/motion_planning/param/costmap_common_params.yaml b/motion_planning/2d/param/costmap_common_params.yaml
similarity index 100%
rename from motion_planning/param/costmap_common_params.yaml
rename to motion_planning/2d/param/costmap_common_params.yaml
diff --git a/motion_planning/param/dwa_local_planner_params.yaml b/motion_planning/2d/param/dwa_local_planner_params.yaml
similarity index 100%
rename from motion_planning/param/dwa_local_planner_params.yaml
rename to motion_planning/2d/param/dwa_local_planner_params.yaml
diff --git a/motion_planning/param/global_costmap_params.yaml b/motion_planning/2d/param/global_costmap_params.yaml
similarity index 100%
rename from motion_planning/param/global_costmap_params.yaml
rename to motion_planning/2d/param/global_costmap_params.yaml
diff --git a/motion_planning/param/global_planner_params.yaml b/motion_planning/2d/param/global_planner_params.yaml
similarity index 100%
rename from motion_planning/param/global_planner_params.yaml
rename to motion_planning/2d/param/global_planner_params.yaml
diff --git a/motion_planning/param/local_costmap_params.yaml b/motion_planning/2d/param/local_costmap_params.yaml
similarity index 100%
rename from motion_planning/param/local_costmap_params.yaml
rename to motion_planning/2d/param/local_costmap_params.yaml
diff --git a/motion_planning/param/move_base_params.yaml b/motion_planning/2d/param/move_base_params.yaml
similarity index 100%
rename from motion_planning/param/move_base_params.yaml
rename to motion_planning/2d/param/move_base_params.yaml
diff --git a/motion_planning/rviz/2d_motion_planning.rviz b/motion_planning/2d/rviz/2d_motion_planning.rviz
similarity index 100%
rename from motion_planning/rviz/2d_motion_planning.rviz
rename to motion_planning/2d/rviz/2d_motion_planning.rviz
diff --git a/motion_planning/3d/ego_planner/ego_rviz.rviz b/motion_planning/3d/ego_planner/ego_rviz.rviz
new file mode 100644
index 0000000..a6d242c
--- /dev/null
+++ b/motion_planning/3d/ego_planner/ego_rviz.rviz
@@ -0,0 +1,188 @@
+Panels:
+ - Class: rviz/Displays
+ Help Height: 78
+ Name: Displays
+ Property Tree Widget:
+ Expanded:
+ - /Global Options1
+ - /Status1
+ Splitter Ratio: 0.5
+ Tree Height: 549
+ - Class: rviz/Selection
+ Name: Selection
+ - Class: rviz/Tool Properties
+ Expanded:
+ - /2D Pose Estimate1
+ - /2D Nav Goal1
+ - /Publish Point1
+ Name: Tool Properties
+ Splitter Ratio: 0.5886790156364441
+ - Class: rviz/Views
+ Expanded:
+ - /Current View1
+ Name: Views
+ Splitter Ratio: 0.5
+ - Class: rviz/Time
+ Experimental: false
+ Name: Time
+ SyncMode: 0
+ SyncSource: PointCloud2
+Preferences:
+ PromptSaveOnExit: true
+Toolbars:
+ toolButtonStyle: 2
+Visualization Manager:
+ Class: ""
+ Displays:
+ - Alpha: 0.5
+ Cell Size: 1
+ Class: rviz/Grid
+ Color: 160; 160; 164
+ Enabled: true
+ Line Style:
+ Line Width: 0.029999999329447746
+ Value: Lines
+ Name: Grid
+ Normal Cell Count: 0
+ Offset:
+ X: 0
+ Y: 0
+ Z: 0
+ Plane: XY
+ Plane Cell Count: 10
+ Reference Frame:
+ Value: true
+ - Class: rviz/Marker
+ Enabled: true
+ Marker Topic: /iris_0_ego_planner_node/optimal_list
+ Name: Marker
+ Namespaces:
+ {}
+ Queue Size: 100
+ Value: true
+ - Alpha: 1
+ Autocompute Intensity Bounds: true
+ Autocompute Value Bounds:
+ Max Value: 10
+ Min Value: -10
+ Value: true
+ Axis: Z
+ Channel Name: intensity
+ Class: rviz/PointCloud2
+ Color: 255; 255; 255
+ Color Transformer: Intensity
+ Decay Time: 0
+ Enabled: true
+ Invert Rainbow: false
+ Max Color: 255; 255; 255
+ Max Intensity: 4096
+ Min Color: 0; 0; 0
+ Min Intensity: 0
+ Name: PointCloud2
+ Position Transformer: XYZ
+ Queue Size: 10
+ Selectable: true
+ Size (Pixels): 3
+ Size (m): 0.009999999776482582
+ Style: Flat Squares
+ Topic: /iris_0_ego_planner_node/grid_map/occupancy_inflate
+ Unreliable: false
+ Use Fixed Frame: true
+ Use rainbow: true
+ Value: true
+ - Alpha: 1
+ Axes Length: 1
+ Axes Radius: 0.10000000149011612
+ Class: rviz/Pose
+ Color: 255; 25; 0
+ Enabled: true
+ Head Length: 0.30000001192092896
+ Head Radius: 0.10000000149011612
+ Name: Pose
+ Shaft Length: 1
+ Shaft Radius: 0.05000000074505806
+ Shape: Arrow
+ Topic: /iris_0/mavros/vision_pose/pose
+ Unreliable: false
+ Value: true
+ - Alpha: 1
+ Axes Length: 1
+ Axes Radius: 0.10000000149011612
+ Class: rviz/Pose
+ Color: 255; 25; 0
+ Enabled: true
+ Head Length: 0.30000001192092896
+ Head Radius: 0.10000000149011612
+ Name: Pose
+ Shaft Length: 1
+ Shaft Radius: 0.05000000074505806
+ Shape: Arrow
+ Topic: /move_base_simple/goal
+ Unreliable: false
+ Value: true
+ Enabled: true
+ Global Options:
+ Background Color: 48; 48; 48
+ Default Light: true
+ Fixed Frame: map
+ Frame Rate: 30
+ Name: root
+ Tools:
+ - Class: rviz/Interact
+ Hide Inactive Objects: true
+ - Class: rviz/MoveCamera
+ - Class: rviz/Select
+ - Class: rviz/FocusCamera
+ - Class: rviz/Measure
+ - Class: rviz/SetInitialPose
+ Theta std deviation: 0.2617993950843811
+ Topic: /initialpose
+ X std deviation: 0.5
+ Y std deviation: 0.5
+ - Class: rviz/SetGoal
+ Topic: /move_base_simple/goal
+ - Class: rviz/PublishPoint
+ Single click: true
+ Topic: /clicked_point
+ Value: true
+ Views:
+ Current:
+ Class: rviz/Orbit
+ Distance: 4.135247707366943
+ Enable Stereo Rendering:
+ Stereo Eye Separation: 0.05999999865889549
+ Stereo Focal Distance: 1
+ Swap Stereo Eyes: false
+ Value: false
+ Focal Point:
+ X: 0.7557446956634521
+ Y: 7.2424845695495605
+ Z: 0.9154675006866455
+ Focal Shape Fixed Size: true
+ Focal Shape Size: 0.05000000074505806
+ Invert Z Axis: false
+ Name: Current View
+ Near Clip Distance: 0.009999999776482582
+ Pitch: 0.8753979206085205
+ Target Frame:
+ Value: Orbit (rviz)
+ Yaw: 3.0754048824310303
+ Saved: ~
+Window Geometry:
+ Displays:
+ collapsed: false
+ Height: 846
+ Hide Left Dock: false
+ Hide Right Dock: false
+ QMainWindow State: 000000ff00000000fd000000040000000000000156000002b0fc0200000008fb0000001200530065006c0065006300740069006f006e00000001e10000009b0000005c00fffffffb0000001e0054006f006f006c002000500072006f007000650072007400690065007302000001ed000001df00000185000000a3fb000000120056006900650077007300200054006f006f02000001df000002110000018500000122fb000000200054006f006f006c002000500072006f0070006500720074006900650073003203000002880000011d000002210000017afb000000100044006900730070006c006100790073010000003d000002b0000000c900fffffffb0000002000730065006c0065006300740069006f006e00200062007500660066006500720200000138000000aa0000023a00000294fb00000014005700690064006500530074006500720065006f02000000e6000000d2000003ee0000030bfb0000000c004b0069006e0065006300740200000186000001060000030c00000261000000010000010f000002b0fc0200000003fb0000001e0054006f006f006c002000500072006f00700065007200740069006500730100000041000000780000000000000000fb0000000a00560069006500770073010000003d000002b0000000a400fffffffb0000001200530065006c0065006300740069006f006e010000025a000000b200000000000000000000000200000490000000a9fc0100000001fb0000000a00560069006500770073030000004e00000080000002e10000019700000003000004b00000003efc0100000002fb0000000800540069006d00650100000000000004b0000002eb00fffffffb0000000800540069006d006501000000000000045000000000000000000000023f000002b000000004000000040000000800000008fc0000000100000002000000010000000a0054006f006f006c00730100000000ffffffff0000000000000000
+ Selection:
+ collapsed: false
+ Time:
+ collapsed: false
+ Tool Properties:
+ collapsed: false
+ Views:
+ collapsed: false
+ Width: 1200
+ X: 67
+ Y: 30
diff --git a/motion_planning/3d/ego_planner/ego_transfer.py b/motion_planning/3d/ego_planner/ego_transfer.py
new file mode 100755
index 0000000..9f1f753
--- /dev/null
+++ b/motion_planning/3d/ego_planner/ego_transfer.py
@@ -0,0 +1,34 @@
+import rospy
+from geometry_msgs.msg import PoseStamped
+import math
+from pyquaternion import Quaternion
+import tf
+import sys
+
+vehicle_type = sys.argv[1]
+local_pose = PoseStamped()
+local_pose.header.frame_id = 'world'
+quaternion = tf.transformations.quaternion_from_euler(-math.pi/2, 0, -math.pi/2)
+q = Quaternion([quaternion[3],quaternion[0],quaternion[1],quaternion[2]])
+
+def vision_callback(data):
+ local_pose.pose.position.x = data.pose.position.x
+ local_pose.pose.position.y = data.pose.position.y
+ local_pose.pose.position.z = data.pose.position.z
+ q_= Quaternion([data.pose.orientation.w,data.pose.orientation.x,data.pose.orientation.y,data.pose.orientation.z])
+ q_ = q_*q
+ local_pose.pose.orientation.w = q_[0]
+ local_pose.pose.orientation.x = q_[1]
+ local_pose.pose.orientation.y = q_[2]
+ local_pose.pose.orientation.z = q_[3]
+
+rospy.init_node('ego_transfer')
+rospy.Subscriber(vehicle_type+"_0/mavros/vision_pose/pose", PoseStamped, vision_callback)
+position_pub = rospy.Publisher(vehicle_type+"_0/camera_pose", PoseStamped, queue_size=2)
+rate = rospy.Rate(60)
+
+while True:
+ local_pose.header.stamp = rospy.Time.now()
+ position_pub.publish(local_pose)
+ rate.sleep()
+
diff --git a/motion_planning/3d/ego_planner/planner/bspline_opt/CMakeLists.txt b/motion_planning/3d/ego_planner/planner/bspline_opt/CMakeLists.txt
new file mode 100755
index 0000000..e61c145
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/bspline_opt/CMakeLists.txt
@@ -0,0 +1,45 @@
+cmake_minimum_required(VERSION 2.8.3)
+project(bspline_opt)
+
+set(CMAKE_BUILD_TYPE "Release")
+ADD_COMPILE_OPTIONS(-std=c++11 )
+ADD_COMPILE_OPTIONS(-std=c++14 )
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -Wall -g")
+
+find_package(Eigen3 REQUIRED)
+find_package(PCL 1.7 REQUIRED)
+
+find_package(catkin REQUIRED COMPONENTS
+ roscpp
+ rospy
+ std_msgs
+ visualization_msgs
+ plan_env
+ cv_bridge
+ traj_utils
+ path_searching
+)
+
+catkin_package(
+ INCLUDE_DIRS include
+ LIBRARIES bspline_opt
+ CATKIN_DEPENDS plan_env
+# DEPENDS system_lib
+)
+
+include_directories(
+ SYSTEM
+ include
+ ${catkin_INCLUDE_DIRS}
+ ${Eigen3_INCLUDE_DIRS}
+ ${PCL_INCLUDE_DIRS}
+)
+
+add_library( bspline_opt
+ src/uniform_bspline.cpp
+ src/bspline_optimizer.cpp
+ src/gradient_descent_optimizer.cpp
+ )
+target_link_libraries( bspline_opt
+ ${catkin_LIBRARIES}
+ )
diff --git a/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/bspline_optimizer.h b/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/bspline_optimizer.h
new file mode 100644
index 0000000..32abae5
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/bspline_optimizer.h
@@ -0,0 +1,212 @@
+#ifndef _BSPLINE_OPTIMIZER_H_
+#define _BSPLINE_OPTIMIZER_H_
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include "bspline_opt/lbfgs.hpp"
+#include
+
+// Gradient and elasitc band optimization
+
+// Input: a signed distance field and a sequence of points
+// Output: the optimized sequence of points
+// The format of points: N x 3 matrix, each row is a point
+namespace ego_planner
+{
+
+ class ControlPoints
+ {
+ public:
+ double clearance;
+ int size;
+ Eigen::MatrixXd points;
+ std::vector> base_point; // The point at the statrt of the direction vector (collision point)
+ std::vector> direction; // Direction vector, must be normalized.
+ std::vector flag_temp; // A flag that used in many places. Initialize it everytime before using it.
+ // std::vector occupancy;
+
+ void resize(const int size_set)
+ {
+ size = size_set;
+
+ base_point.clear();
+ direction.clear();
+ flag_temp.clear();
+ // occupancy.clear();
+
+ points.resize(3, size_set);
+ base_point.resize(size);
+ direction.resize(size);
+ flag_temp.resize(size);
+ // occupancy.resize(size);
+ }
+
+ void segment(ControlPoints &buf, const int start, const int end)
+ {
+ if (start < 0 || end >= size || points.rows() != 3)
+ {
+ ROS_ERROR("Wrong segment index! start=%d, end=%d", start, end);
+ return;
+ }
+
+ buf.resize(end - start + 1);
+ buf.points = points.block(0, start, 3, end - start + 1);
+ buf.clearance = clearance;
+ buf.size = end - start + 1;
+ for (int i = start; i <= end; i++)
+ {
+ buf.base_point[i - start] = base_point[i];
+ buf.direction[i - start] = direction[i];
+
+ // if ( buf.base_point[i - start].size() > 1 )
+ // {
+ // ROS_ERROR("buf.base_point[i - start].size()=%d, base_point[i].size()=%d", buf.base_point[i - start].size(), base_point[i].size());
+ // }
+ }
+
+ // cout << "RichInfoOneSeg_temp, insede" << endl;
+ // for ( int k=0; k 0 )
+ // {
+ // cout << "###" << buf.points.col(k).transpose() << endl;
+ // for (int k2 = 0; k2 < buf.base_point[k].size(); k2++)
+ // {
+ // cout << " " << buf.base_point[k][k2].transpose() << " @ " << buf.direction[k][k2].transpose() << endl;
+ // }
+ // }
+ }
+ };
+
+ class BsplineOptimizer
+ {
+
+ public:
+ BsplineOptimizer() {}
+ ~BsplineOptimizer() {}
+
+ /* main API */
+ void setEnvironment(const GridMap::Ptr &map);
+ void setEnvironment(const GridMap::Ptr &map, const fast_planner::ObjPredictor::Ptr mov_obj);
+ void setParam(ros::NodeHandle &nh);
+ Eigen::MatrixXd BsplineOptimizeTraj(const Eigen::MatrixXd &points, const double &ts,
+ const int &cost_function, int max_num_id, int max_time_id);
+
+ /* helper function */
+
+ // required inputs
+ void setControlPoints(const Eigen::MatrixXd &points);
+ void setBsplineInterval(const double &ts);
+ void setSwarmTrajs(SwarmTrajData *swarm_trajs_ptr);
+ void setDroneId(const int drone_id);
+
+ // optional inputs
+ void setGuidePath(const vector &guide_pt);
+ void setWaypoints(const vector &waypts,
+ const vector &waypt_idx); // N-2 constraints at most
+ void setLocalTargetPt(const Eigen::Vector3d local_target_pt) { local_target_pt_ = local_target_pt; };
+
+ void optimize();
+
+ ControlPoints getControlPoints() { return cps_; };
+
+ AStar::Ptr a_star_;
+ std::vector ref_pts_;
+
+ std::vector distinctiveTrajs(vector> segments);
+ std::vector> initControlPoints(Eigen::MatrixXd &init_points, bool flag_first_init = true);
+ bool BsplineOptimizeTrajRebound(Eigen::MatrixXd &optimal_points, double ts); // must be called after initControlPoints()
+ bool BsplineOptimizeTrajRebound(Eigen::MatrixXd &optimal_points, double &final_cost, const ControlPoints &control_points, double ts);
+ bool BsplineOptimizeTrajRefine(const Eigen::MatrixXd &init_points, const double ts, Eigen::MatrixXd &optimal_points);
+
+ inline int getOrder(void) { return order_; }
+ inline double getSwarmClearance(void) { return swarm_clearance_; }
+
+ private:
+ GridMap::Ptr grid_map_;
+ fast_planner::ObjPredictor::Ptr moving_objs_;
+ SwarmTrajData *swarm_trajs_{NULL}; // Can not use shared_ptr and no need to free
+ int drone_id_;
+
+ enum FORCE_STOP_OPTIMIZE_TYPE
+ {
+ DONT_STOP,
+ STOP_FOR_REBOUND,
+ STOP_FOR_ERROR
+ } force_stop_type_;
+
+ // main input
+ // Eigen::MatrixXd control_points_; // B-spline control points, N x dim
+ double bspline_interval_; // B-spline knot span
+ Eigen::Vector3d end_pt_; // end of the trajectory
+ // int dim_; // dimension of the B-spline
+ //
+ vector guide_pts_; // geometric guiding path points, N-6
+ vector waypoints_; // waypts constraints
+ vector waypt_idx_; // waypts constraints index
+ //
+ int max_num_id_, max_time_id_; // stopping criteria
+ int cost_function_; // used to determine objective function
+ double start_time_; // global time for moving obstacles
+
+ /* optimization parameters */
+ int order_; // bspline degree
+ double lambda1_; // jerk smoothness weight
+ double lambda2_, new_lambda2_; // distance weight
+ double lambda3_; // feasibility weight
+ double lambda4_; // curve fitting
+
+ int a;
+ //
+ double dist0_, swarm_clearance_; // safe distance
+ double max_vel_, max_acc_; // dynamic limits
+
+ int variable_num_; // optimization variables
+ int iter_num_; // iteration of the solver
+ Eigen::VectorXd best_variable_; //
+ double min_cost_; //
+
+ Eigen::Vector3d local_target_pt_;
+
+#define INIT_min_ellip_dist_ 123456789.0123456789
+ double min_ellip_dist_;
+
+ ControlPoints cps_;
+
+ /* cost function */
+ /* calculate each part of cost function with control points q as input */
+
+ static double costFunction(const std::vector &x, std::vector &grad, void *func_data);
+ void combineCost(const std::vector &x, vector &grad, double &cost);
+
+ // q contains all control points
+ void calcSmoothnessCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient, bool falg_use_jerk = true);
+ void calcFeasibilityCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient);
+ void calcTerminalCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient);
+ void calcDistanceCostRebound(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient, int iter_num, double smoothness_cost);
+ void calcMovingObjCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient);
+ void calcSwarmCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient);
+ void calcFitnessCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient);
+ bool check_collision_and_rebound(void);
+
+ static int earlyExit(void *func_data, const double *x, const double *g, const double fx, const double xnorm, const double gnorm, const double step, int n, int k, int ls);
+ static double costFunctionRebound(void *func_data, const double *x, double *grad, const int n);
+ static double costFunctionRefine(void *func_data, const double *x, double *grad, const int n);
+
+ bool rebound_optimize(double &final_cost);
+ bool refine_optimize();
+ void combineCostRebound(const double *x, double *grad, double &f_combine, const int n);
+ void combineCostRefine(const double *x, double *grad, double &f_combine, const int n);
+
+ /* for benckmark evaluation only */
+ public:
+ typedef unique_ptr Ptr;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+ };
+
+} // namespace ego_planner
+#endif
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/gradient_descent_optimizer.h b/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/gradient_descent_optimizer.h
new file mode 100755
index 0000000..01a5a71
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/gradient_descent_optimizer.h
@@ -0,0 +1,52 @@
+#ifndef _GRADIENT_DESCENT_OPT_H_
+#define _GRADIENT_DESCENT_OPT_H_
+
+#include
+#include
+#include
+
+using namespace std;
+
+class GradientDescentOptimizer
+{
+
+public:
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+
+ typedef double (*objfunDef)(const Eigen::VectorXd &x, Eigen::VectorXd &grad, bool &force_return, void *data);
+ enum RESULT
+ {
+ FIND_MIN,
+ FAILED,
+ RETURN_BY_ORDER,
+ REACH_MAX_ITERATION
+ };
+
+ GradientDescentOptimizer(int v_num, objfunDef objf, void *f_data)
+ {
+ variable_num_ = v_num;
+ objfun_ = objf;
+ f_data_ = f_data;
+ };
+
+ void set_maxiter(int limit) { iter_limit_ = limit; }
+ void set_maxeval(int limit) { invoke_limit_ = limit; }
+ void set_xtol_rel(double xtol_rel) { xtol_rel_ = xtol_rel; }
+ void set_xtol_abs(double xtol_abs) { xtol_abs_ = xtol_abs; }
+ void set_min_grad(double min_grad) { min_grad_ = min_grad; }
+
+ RESULT optimize(Eigen::VectorXd &x_init_optimal, double &opt_f);
+
+private:
+ int variable_num_{0};
+ int iter_limit_{1e10};
+ int invoke_limit_{1e10};
+ double xtol_rel_;
+ double xtol_abs_;
+ double min_grad_;
+ double time_limit_;
+ void *f_data_;
+ objfunDef objfun_;
+};
+
+#endif
diff --git a/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/lbfgs.hpp b/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/lbfgs.hpp
new file mode 100644
index 0000000..484c4b6
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/lbfgs.hpp
@@ -0,0 +1,1451 @@
+#ifndef LBFGS_HPP
+#define LBFGS_HPP
+
+#include
+#include
+#include
+#include
+#include
+
+namespace lbfgs
+{
+ // ----------------------- Data Type Part -----------------------
+
+ /**
+ * Return values of lbfgs_optimize().
+ *
+ * Roughly speaking, a negative value indicates an error.
+ */
+
+ enum
+ {
+ /** L-BFGS reaches convergence. */
+ LBFGS_CONVERGENCE = 0,
+ /** L-BFGS satisfies stopping criteria. */
+ LBFGS_STOP,
+ /** The initial variables already minimize the objective function. */
+ LBFGS_ALREADY_MINIMIZED,
+
+ /** Unknown error. */
+ LBFGSERR_UNKNOWNERROR = -1024,
+ /** Logic error. */
+ LBFGSERR_LOGICERROR,
+ /** The minimization process has been canceled. */
+ LBFGSERR_CANCELED,
+ /** Invalid number of variables specified. */
+ LBFGSERR_INVALID_N,
+ /** Invalid parameter lbfgs_parameter_t::mem_size specified. */
+ LBFGSERR_INVALID_MEMSIZE,
+ /** Invalid parameter lbfgs_parameter_t::g_epsilon specified. */
+ LBFGSERR_INVALID_GEPSILON,
+ /** Invalid parameter lbfgs_parameter_t::past specified. */
+ LBFGSERR_INVALID_TESTPERIOD,
+ /** Invalid parameter lbfgs_parameter_t::delta specified. */
+ LBFGSERR_INVALID_DELTA,
+ /** Invalid parameter lbfgs_parameter_t::min_step specified. */
+ LBFGSERR_INVALID_MINSTEP,
+ /** Invalid parameter lbfgs_parameter_t::max_step specified. */
+ LBFGSERR_INVALID_MAXSTEP,
+ /** Invalid parameter lbfgs_parameter_t::f_dec_coeff specified. */
+ LBFGSERR_INVALID_FDECCOEFF,
+ /** Invalid parameter lbfgs_parameter_t::wolfe specified. */
+ LBFGSERR_INVALID_SCURVCOEFF,
+ /** Invalid parameter lbfgs_parameter_t::xtol specified. */
+ LBFGSERR_INVALID_XTOL,
+ /** Invalid parameter lbfgs_parameter_t::max_linesearch specified. */
+ LBFGSERR_INVALID_MAXLINESEARCH,
+ /** The line-search step went out of the interval of uncertainty. */
+ LBFGSERR_OUTOFINTERVAL,
+ /** A logic error occurred; alternatively, the interval of uncertainty
+ became too small. */
+ LBFGSERR_INCORRECT_TMINMAX,
+ /** A rounding error occurred; alternatively, no line-search step
+ satisfies the sufficient decrease and curvature conditions. */
+ LBFGSERR_ROUNDING_ERROR,
+ /** The line-search step became smaller than lbfgs_parameter_t::min_step. */
+ LBFGSERR_MINIMUMSTEP,
+ /** The line-search step became larger than lbfgs_parameter_t::max_step. */
+ LBFGSERR_MAXIMUMSTEP,
+ /** The line-search routine reaches the maximum number of evaluations. */
+ LBFGSERR_MAXIMUMLINESEARCH,
+ /** The algorithm routine reaches the maximum number of iterations. */
+ LBFGSERR_MAXIMUMITERATION,
+ /** Relative width of the interval of uncertainty is at most
+ lbfgs_parameter_t::xtol. */
+ LBFGSERR_WIDTHTOOSMALL,
+ /** A logic error (negative line-search step) occurred. */
+ LBFGSERR_INVALIDPARAMETERS,
+ /** The current search direction increases the objective function value. */
+ LBFGSERR_INCREASEGRADIENT,
+ };
+
+ /**
+ * L-BFGS optimization parameters.
+ * Call lbfgs_load_default_parameters() function to initialize parameters to the
+ * default values.
+ */
+ struct lbfgs_parameter_t
+ {
+ /**
+ * The number of corrections to approximate the inverse hessian matrix.
+ * The L-BFGS routine stores the computation results of previous m
+ * iterations to approximate the inverse hessian matrix of the current
+ * iteration. This parameter controls the size of the limited memories
+ * (corrections). The default value is 6. Values less than 3 are
+ * not recommended. Large values will result in excessive computing time.
+ */
+ int mem_size;
+
+ /**
+ * Epsilon for grad norm convergence test.
+ * This parameter determines the accuracy with which the solution is to
+ * be found. A minimization terminates when
+ * ||g|| < g_epsilon * max(1, ||x||),
+ * where ||.|| denotes the Euclidean (L2) norm. The default value is 1e-5.
+ */
+ double g_epsilon;
+
+ /**
+ * Distance for delta-based convergence test.
+ * This parameter determines the distance, in iterations, to compute
+ * the rate of decrease of the objective function. If the value of this
+ * parameter is zero, the library does not perform the delta-based
+ * convergence test. The default value is 0.
+ */
+ int past;
+
+ /**
+ * Delta for convergence test.
+ * This parameter determines the minimum rate of decrease of the
+ * objective function. The library stops iterations when the
+ * following condition is met:
+ * (f' - f) / f < delta,
+ * where f' is the objective value of past iterations ago, and f is
+ * the objective value of the current iteration.
+ * The default value is 1e-5.
+ */
+ double delta;
+
+ /**
+ * The maximum number of iterations.
+ * The lbfgs_optimize() function terminates an optimization process with
+ * ::LBFGSERR_MAXIMUMITERATION status code when the iteration count
+ * exceedes this parameter. Setting this parameter to zero continues an
+ * optimization process until a convergence or error. The default value
+ * is 0.
+ */
+ int max_iterations;
+
+ /**
+ * The maximum number of trials for the line search.
+ * This parameter controls the number of function and gradients evaluations
+ * per iteration for the line search routine. The default value is 40.
+ */
+ int max_linesearch;
+
+ /**
+ * The minimum step of the line search routine.
+ * The default value is 1e-20. This value need not be modified unless
+ * the exponents are too large for the machine being used, or unless the
+ * problem is extremely badly scaled (in which case the exponents should
+ * be increased).
+ */
+ double min_step;
+
+ /**
+ * The maximum step of the line search.
+ * The default value is 1e+20. This value need not be modified unless
+ * the exponents are too large for the machine being used, or unless the
+ * problem is extremely badly scaled (in which case the exponents should
+ * be increased).
+ */
+ double max_step;
+
+ /**
+ * A parameter to control the accuracy of the line search routine.
+ * The default value is 1e-4. This parameter should be greater
+ * than zero and smaller than 0.5.
+ */
+ double f_dec_coeff;
+
+ /**
+ * A parameter to control the accuracy of the line search routine.
+ * The default value is 0.9. If the function and gradient
+ * evaluations are inexpensive with respect to the cost of the
+ * iteration (which is sometimes the case when solving very large
+ * problems) it may be advantageous to set this parameter to a small
+ * value. A typical small value is 0.1. This parameter shuold be
+ * greater than the f_dec_coeff parameter (1e-4)
+ * and smaller than 1.0.
+ */
+ double s_curv_coeff;
+
+ /**
+ * The machine precision for floating-point values.
+ * This parameter must be a positive value set by a client program to
+ * estimate the machine precision. The line search routine will terminate
+ * with the status code (::LBFGSERR_ROUNDING_ERROR) if the relative width
+ * of the interval of uncertainty is less than this parameter.
+ */
+ double xtol;
+ };
+
+ /**
+ * Callback interface to provide objective function and gradient evaluations.
+ *
+ * The lbfgs_optimize() function call this function to obtain the values of objective
+ * function and its gradients when needed. A client program must implement
+ * this function to evaluate the values of the objective function and its
+ * gradients, given current values of variables.
+ *
+ * @param instance The user data sent for lbfgs_optimize() function by the client.
+ * @param x The current values of variables.
+ * @param g The gradient vector. The callback function must compute
+ * the gradient values for the current variables.
+ * @param n The number of variables.
+ * @retval double The value of the objective function for the current
+ * variables.
+ */
+ typedef double (*lbfgs_evaluate_t)(void *instance,
+ const double *x,
+ double *g,
+ const int n);
+
+ /**
+ * Callback interface to provide an upper bound at the beginning of current linear search.
+ *
+ * The lbfgs_optimize() function call this function to obtain the values of the
+ * upperbound of the stepsize to search in, provided with the beginning values of
+ * variables before the linear search, and the current step vector (can be descent direction).
+ * A client program can implement this function for more efficient linesearch.
+ * If it is not used, just set it NULL or nullptr.
+ *
+ * @param instance The user data sent for lbfgs_optimize() function by the client.
+ * @param xp The values of variables before current line search.
+ * @param d The step vector. It can be the descent direction.
+ * @param n The number of variables.
+ * @retval double The upperboud of the step in current line search routine,
+ * such that stpbound*d is the maximum reasonable step.
+ */
+ typedef double (*lbfgs_stepbound_t)(void *instance,
+ const double *xp,
+ const double *d,
+ const int n);
+
+ /**
+ * Callback interface to receive the progress of the optimization process.
+ *
+ * The lbfgs_optimize() function call this function for each iteration. Implementing
+ * this function, a client program can store or display the current progress
+ * of the optimization process. If it is not used, just set it NULL or nullptr.
+ *
+ * @param instance The user data sent for lbfgs_optimize() function by the client.
+ * @param x The current values of variables.
+ * @param g The current gradient values of variables.
+ * @param fx The current value of the objective function.
+ * @param xnorm The Euclidean norm of the variables.
+ * @param gnorm The Euclidean norm of the gradients.
+ * @param step The line-search step used for this iteration.
+ * @param n The number of variables.
+ * @param k The iteration count.
+ * @param ls The number of evaluations called for this iteration.
+ * @retval int Zero to continue the optimization process. Returning a
+ * non-zero value will cancel the optimization process.
+ */
+ typedef int (*lbfgs_progress_t)(void *instance,
+ const double *x,
+ const double *g,
+ const double fx,
+ const double xnorm,
+ const double gnorm,
+ const double step,
+ int n,
+ int k,
+ int ls);
+
+ /**
+ * Callback data struct
+ */
+
+ struct callback_data_t
+ {
+ int n;
+ void *instance;
+ lbfgs_evaluate_t proc_evaluate;
+ lbfgs_stepbound_t proc_stepbound;
+ lbfgs_progress_t proc_progress;
+ };
+
+ /**
+ * Iteration data struct
+ */
+ struct iteration_data_t
+ {
+ double alpha;
+ double *s; /* [n] */
+ double *y; /* [n] */
+ double ys; /* vecdot(y, s) */
+ };
+
+ // ----------------------- Arithmetic Part -----------------------
+
+/**
+ * Define the local variables for computing minimizers.
+ */
+#define USES_MINIMIZER_LBFGS \
+ double a, d, gamm, theta, p, q, r, s;
+
+/**
+ * Find a minimizer of an interpolated cubic function.
+ * @param cm The minimizer of the interpolated cubic.
+ * @param u The value of one point, u.
+ * @param fu The value of f(u).
+ * @param du The value of f'(u).
+ * @param v The value of another point, v.
+ * @param fv The value of f(v).
+ * @param du The value of f'(v).
+ */
+#define CUBIC_MINIMIZER_LBFGS(cm, u, fu, du, v, fv, dv) \
+ d = (v) - (u); \
+ theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \
+ p = fabs(theta); \
+ q = fabs(du); \
+ r = fabs(dv); \
+ s = p >= q ? p : q; \
+ s = s >= r ? s : r; \
+ /* gamm = s*sqrt((theta/s)**2 - (du/s) * (dv/s)) */ \
+ a = theta / s; \
+ gamm = s * sqrt(a * a - ((du) / s) * ((dv) / s)); \
+ if ((v) < (u)) \
+ gamm = -gamm; \
+ p = gamm - (du) + theta; \
+ q = gamm - (du) + gamm + (dv); \
+ r = p / q; \
+ (cm) = (u) + r * d;
+
+/**
+ * Find a minimizer of an interpolated cubic function.
+ * @param cm The minimizer of the interpolated cubic.
+ * @param u The value of one point, u.
+ * @param fu The value of f(u).
+ * @param du The value of f'(u).
+ * @param v The value of another point, v.
+ * @param fv The value of f(v).
+ * @param du The value of f'(v).
+ * @param xmin The maximum value.
+ * @param xmin The minimum value.
+ */
+#define CUBIC_MINIMIZER2_LBFGS(cm, u, fu, du, v, fv, dv, xmin, xmax) \
+ d = (v) - (u); \
+ theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \
+ p = fabs(theta); \
+ q = fabs(du); \
+ r = fabs(dv); \
+ s = p >= q ? p : q; \
+ s = s >= r ? s : r; \
+ /* gamm = s*sqrt((theta/s)**2 - (du/s) * (dv/s)) */ \
+ a = theta / s; \
+ gamm = a * a - ((du) / s) * ((dv) / s); \
+ gamm = gamm > 0 ? s * sqrt(gamm) : 0; \
+ if ((u) < (v)) \
+ gamm = -gamm; \
+ p = gamm - (dv) + theta; \
+ q = gamm - (dv) + gamm + (du); \
+ r = p / q; \
+ if (r < 0. && gamm != 0.) \
+ { \
+ (cm) = (v)-r * d; \
+ } \
+ else if (a < 0) \
+ { \
+ (cm) = (xmax); \
+ } \
+ else \
+ { \
+ (cm) = (xmin); \
+ }
+
+/**
+ * Find a minimizer of an interpolated quadratic function.
+ * @param qm The minimizer of the interpolated quadratic.
+ * @param u The value of one point, u.
+ * @param fu The value of f(u).
+ * @param du The value of f'(u).
+ * @param v The value of another point, v.
+ * @param fv The value of f(v).
+ */
+#define QUARD_MINIMIZER_LBFGS(qm, u, fu, du, v, fv) \
+ a = (v) - (u); \
+ (qm) = (u) + (du) / (((fu) - (fv)) / a + (du)) / 2 * a;
+
+/**
+ * Find a minimizer of an interpolated quadratic function.
+ * @param qm The minimizer of the interpolated quadratic.
+ * @param u The value of one point, u.
+ * @param du The value of f'(u).
+ * @param v The value of another point, v.
+ * @param dv The value of f'(v).
+ */
+#define QUARD_MINIMIZER2_LBFGS(qm, u, du, v, dv) \
+ a = (u) - (v); \
+ (qm) = (v) + (dv) / ((dv) - (du)) * a;
+
+ inline void *vecalloc(size_t size)
+ {
+ void *memblock = malloc(size);
+ if (memblock)
+ {
+ memset(memblock, 0, size);
+ }
+ return memblock;
+ }
+
+ inline void vecfree(void *memblock)
+ {
+ free(memblock);
+ }
+
+ inline void veccpy(double *y, const double *x, const int n)
+ {
+ memcpy(y, x, sizeof(double) * n);
+ }
+
+ inline void vecncpy(double *y, const double *x, const int n)
+ {
+ int i;
+
+ for (i = 0; i < n; ++i)
+ {
+ y[i] = -x[i];
+ }
+ }
+
+ inline void vecadd(double *y, const double *x, const double c, const int n)
+ {
+ int i;
+
+ for (i = 0; i < n; ++i)
+ {
+ y[i] += c * x[i];
+ }
+ }
+
+ inline void vecdiff(double *z, const double *x, const double *y, const int n)
+ {
+ int i;
+
+ for (i = 0; i < n; ++i)
+ {
+ z[i] = x[i] - y[i];
+ }
+ }
+
+ inline void vecscale(double *y, const double c, const int n)
+ {
+ int i;
+
+ for (i = 0; i < n; ++i)
+ {
+ y[i] *= c;
+ }
+ }
+
+ inline void vecdot(double *s, const double *x, const double *y, const int n)
+ {
+ int i;
+ *s = 0.;
+ for (i = 0; i < n; ++i)
+ {
+ *s += x[i] * y[i];
+ }
+ }
+
+ inline void vec2norm(double *s, const double *x, const int n)
+ {
+ vecdot(s, x, x, n);
+ *s = (double)sqrt(*s);
+ }
+
+ inline void vec2norminv(double *s, const double *x, const int n)
+ {
+ vec2norm(s, x, n);
+ *s = (double)(1.0 / *s);
+ }
+
+ // ----------------------- L-BFGS Part -----------------------
+
+ /**
+ * Update a safeguarded trial value and interval for line search.
+ *
+ * The parameter x represents the step with the least function value.
+ * The parameter t represents the current step. This function assumes
+ * that the derivative at the point of x in the direction of the step.
+ * If the bracket is set to true, the minimizer has been bracketed in
+ * an interval of uncertainty with endpoints between x and y.
+ *
+ * @param x The pointer to the value of one endpoint.
+ * @param fx The pointer to the value of f(x).
+ * @param dx The pointer to the value of f'(x).
+ * @param y The pointer to the value of another endpoint.
+ * @param fy The pointer to the value of f(y).
+ * @param dy The pointer to the value of f'(y).
+ * @param t The pointer to the value of the trial value, t.
+ * @param ft The pointer to the value of f(t).
+ * @param dt The pointer to the value of f'(t).
+ * @param tmin The minimum value for the trial value, t.
+ * @param tmax The maximum value for the trial value, t.
+ * @param brackt The pointer to the predicate if the trial value is
+ * bracketed.
+ * @retval int Status value. Zero indicates a normal termination.
+ *
+ * @see
+ * Jorge J. More and David J. Thuente. Line search algorithm with
+ * guaranteed sufficient decrease. ACM Transactions on Mathematical
+ * Software (TOMS), Vol 20, No 3, pp. 286-307, 1994.
+ */
+ inline int update_trial_interval(double *x,
+ double *fx,
+ double *dx,
+ double *y,
+ double *fy,
+ double *dy,
+ double *t,
+ double *ft,
+ double *dt,
+ const double tmin,
+ const double tmax,
+ int *brackt)
+ {
+ int bound;
+ int dsign = *dt * (*dx / fabs(*dx)) < 0.;
+ double mc; /* minimizer of an interpolated cubic. */
+ double mq; /* minimizer of an interpolated quadratic. */
+ double newt; /* new trial value. */
+ USES_MINIMIZER_LBFGS; /* for CUBIC_MINIMIZER and QUARD_MINIMIZER. */
+
+ /* Check the input parameters for errors. */
+ if (*brackt)
+ {
+ if (*t <= (*x <= *y ? *x : *y) || (*x >= *y ? *x : *y) <= *t)
+ {
+ /* The trival value t is out of the interval. */
+ return LBFGSERR_OUTOFINTERVAL;
+ }
+ if (0. <= *dx * (*t - *x))
+ {
+ /* The function must decrease from x. */
+ return LBFGSERR_INCREASEGRADIENT;
+ }
+ if (tmax < tmin)
+ {
+ /* Incorrect tmin and tmax specified. */
+ return LBFGSERR_INCORRECT_TMINMAX;
+ }
+ }
+
+ /*
+ Trial value selection.
+ */
+ if (*fx < *ft)
+ {
+ /*
+ Case 1: a higher function value.
+ The minimum is brackt. If the cubic minimizer is closer
+ to x than the quadratic one, the cubic one is taken, else
+ the average of the minimizers is taken.
+ */
+ *brackt = 1;
+ bound = 1;
+ CUBIC_MINIMIZER_LBFGS(mc, *x, *fx, *dx, *t, *ft, *dt);
+ QUARD_MINIMIZER_LBFGS(mq, *x, *fx, *dx, *t, *ft);
+ if (fabs(mc - *x) < fabs(mq - *x))
+ {
+ newt = mc;
+ }
+ else
+ {
+ newt = mc + 0.5 * (mq - mc);
+ }
+ }
+ else if (dsign)
+ {
+ /*
+ Case 2: a lower function value and derivatives of
+ opposite sign. The minimum is brackt. If the cubic
+ minimizer is closer to x than the quadratic (secant) one,
+ the cubic one is taken, else the quadratic one is taken.
+ */
+ *brackt = 1;
+ bound = 0;
+ CUBIC_MINIMIZER_LBFGS(mc, *x, *fx, *dx, *t, *ft, *dt);
+ QUARD_MINIMIZER2_LBFGS(mq, *x, *dx, *t, *dt);
+ if (fabs(mc - *t) > fabs(mq - *t))
+ {
+ newt = mc;
+ }
+ else
+ {
+ newt = mq;
+ }
+ }
+ else if (fabs(*dt) < fabs(*dx))
+ {
+ /*
+ Case 3: a lower function value, derivatives of the
+ same sign, and the magnitude of the derivative decreases.
+ The cubic minimizer is only used if the cubic tends to
+ infinity in the direction of the minimizer or if the minimum
+ of the cubic is beyond t. Otherwise the cubic minimizer is
+ defined to be either tmin or tmax. The quadratic (secant)
+ minimizer is also computed and if the minimum is brackt
+ then the the minimizer closest to x is taken, else the one
+ farthest away is taken.
+ */
+ bound = 1;
+ CUBIC_MINIMIZER2_LBFGS(mc, *x, *fx, *dx, *t, *ft, *dt, tmin, tmax);
+ QUARD_MINIMIZER2_LBFGS(mq, *x, *dx, *t, *dt);
+ if (*brackt)
+ {
+ if (fabs(*t - mc) < fabs(*t - mq))
+ {
+ newt = mc;
+ }
+ else
+ {
+ newt = mq;
+ }
+ }
+ else
+ {
+ if (fabs(*t - mc) > fabs(*t - mq))
+ {
+ newt = mc;
+ }
+ else
+ {
+ newt = mq;
+ }
+ }
+ }
+ else
+ {
+ /*
+ Case 4: a lower function value, derivatives of the
+ same sign, and the magnitude of the derivative does
+ not decrease. If the minimum is not brackt, the step
+ is either tmin or tmax, else the cubic minimizer is taken.
+ */
+ bound = 0;
+ if (*brackt)
+ {
+ CUBIC_MINIMIZER_LBFGS(newt, *t, *ft, *dt, *y, *fy, *dy);
+ }
+ else if (*x < *t)
+ {
+ newt = tmax;
+ }
+ else
+ {
+ newt = tmin;
+ }
+ }
+
+ /*
+ Update the interval of uncertainty. This update does not
+ depend on the new step or the case analysis above.
+
+ - Case a: if f(x) < f(t),
+ x <- x, y <- t.
+ - Case b: if f(t) <= f(x) && f'(t)*f'(x) > 0,
+ x <- t, y <- y.
+ - Case c: if f(t) <= f(x) && f'(t)*f'(x) < 0,
+ x <- t, y <- x.
+ */
+ if (*fx < *ft)
+ {
+ /* Case a */
+ *y = *t;
+ *fy = *ft;
+ *dy = *dt;
+ }
+ else
+ {
+ /* Case c */
+ if (dsign)
+ {
+ *y = *x;
+ *fy = *fx;
+ *dy = *dx;
+ }
+ /* Cases b and c */
+ *x = *t;
+ *fx = *ft;
+ *dx = *dt;
+ }
+
+ /* Clip the new trial value in [tmin, tmax]. */
+ if (tmax < newt)
+ newt = tmax;
+ if (newt < tmin)
+ newt = tmin;
+
+ /*
+ Redefine the new trial value if it is close to the upper bound
+ of the interval.
+ */
+ if (*brackt && bound)
+ {
+ mq = *x + 0.66 * (*y - *x);
+ if (*x < *y)
+ {
+ if (mq < newt)
+ newt = mq;
+ }
+ else
+ {
+ if (newt < mq)
+ newt = mq;
+ }
+ }
+
+ /* Return the new trial value. */
+ *t = newt;
+ return 0;
+ }
+
+ inline int line_search_morethuente(int n,
+ double *x,
+ double *f,
+ double *g,
+ double *s,
+ double *stp,
+ const double *xp,
+ const double *gp,
+ const double *stpmin,
+ const double *stpmax,
+ callback_data_t *cd,
+ const lbfgs_parameter_t *param)
+ {
+ int count = 0;
+ int brackt, stage1, uinfo = 0;
+ double dg;
+ double stx, fx, dgx;
+ double sty, fy, dgy;
+ double fxm, dgxm, fym, dgym, fm, dgm;
+ double finit, ftest1, dginit, dgtest;
+ double width, prev_width;
+ double stmin, stmax;
+
+ /* Check the input parameters for errors. */
+ if (*stp <= 0.)
+ {
+ return LBFGSERR_INVALIDPARAMETERS;
+ }
+
+ /* Compute the initial gradient in the search direction. */
+ vecdot(&dginit, g, s, n);
+
+ /* Make sure that s points to a descent direction. */
+ if (0 < dginit)
+ {
+ return LBFGSERR_INCREASEGRADIENT;
+ }
+
+ /* Initialize local variables. */
+ brackt = 0;
+ stage1 = 1;
+ finit = *f;
+ dgtest = param->f_dec_coeff * dginit;
+ width = *stpmax - *stpmin;
+ prev_width = 2.0 * width;
+
+ /*
+ The variables stx, fx, dgx contain the values of the step,
+ function, and directional derivative at the best step.
+ The variables sty, fy, dgy contain the value of the step,
+ function, and derivative at the other endpoint of
+ the interval of uncertainty.
+ The variables stp, f, dg contain the values of the step,
+ function, and derivative at the current step.
+ */
+ stx = sty = 0.;
+ fx = fy = finit;
+ dgx = dgy = dginit;
+
+ for (;;)
+ {
+ /* Report the progress. */
+ if (cd->proc_progress)
+ {
+ double xnorm;
+ double gnorm;
+ vec2norm(&xnorm, x, n);
+ vec2norm(&gnorm, g, n);
+ if (cd->proc_progress(cd->instance, x, g, fx, xnorm, gnorm, *stp, cd->n, 0, 0))
+ {
+ return LBFGSERR_CANCELED;
+ }
+ }
+
+ /*
+ Set the minimum and maximum steps to correspond to the
+ present interval of uncertainty.
+ */
+ if (brackt)
+ {
+ stmin = stx <= sty ? stx : sty;
+ stmax = stx >= sty ? stx : sty;
+ }
+ else
+ {
+ stmin = stx;
+ stmax = *stp + 4.0 * (*stp - stx);
+ }
+
+ /* Clip the step in the range of [stpmin, stpmax]. */
+ if (*stp < *stpmin)
+ *stp = *stpmin;
+ if (*stpmax < *stp)
+ *stp = *stpmax;
+
+ /*
+ If an unusual termination is to occur then let
+ stp be the lowest point obtained so far.
+ */
+ if ((brackt && ((*stp <= stmin || stmax <= *stp) || param->max_linesearch <= count + 1 || uinfo != 0)) || (brackt && (stmax - stmin <= param->xtol * stmax)))
+ {
+ *stp = stx;
+ }
+
+ /*
+ Compute the current value of x:
+ x <- x + (*stp) * s.
+ */
+ veccpy(x, xp, n);
+ vecadd(x, s, *stp, n);
+
+ /* Evaluate the function and gradient values. */
+ *f = cd->proc_evaluate(cd->instance, x, g, cd->n);
+ vecdot(&dg, g, s, n);
+
+ ftest1 = finit + *stp * dgtest;
+ ++count;
+
+ /* Test for errors and convergence. */
+ if (brackt && ((*stp <= stmin || stmax <= *stp) || uinfo != 0))
+ {
+ /* Rounding errors prevent further progress. */
+ return LBFGSERR_ROUNDING_ERROR;
+ }
+ if (*stp == *stpmax && *f <= ftest1 && dg <= dgtest)
+ {
+ /* The step is the maximum value. */
+ return LBFGSERR_MAXIMUMSTEP;
+ }
+ if (*stp == *stpmin && (ftest1 < *f || dgtest <= dg))
+ {
+ /* The step is the minimum value. */
+ return LBFGSERR_MINIMUMSTEP;
+ }
+ if (brackt && (stmax - stmin) <= param->xtol * stmax)
+ {
+ /* Relative width of the interval of uncertainty is at most xtol. */
+ return LBFGSERR_WIDTHTOOSMALL;
+ }
+ if (param->max_linesearch <= count)
+ {
+ /* Maximum number of iteration. */
+ return LBFGSERR_MAXIMUMLINESEARCH;
+ }
+ if (*f <= ftest1 && fabs(dg) <= param->s_curv_coeff * (-dginit))
+ {
+ /* The sufficient decrease condition and the strong curvature condition hold. */
+ return count;
+ }
+
+ /*
+ In the first stage we seek a step for which the modified
+ function has a nonpositive value and nonnegative derivative.
+ */
+ if (stage1 && *f <= ftest1 &&
+ (param->f_dec_coeff <= param->s_curv_coeff ? param->f_dec_coeff : param->s_curv_coeff) * dginit <= dg)
+ {
+ stage1 = 0;
+ }
+
+ /*
+ A modified function is used to predict the step only if
+ we have not obtained a step for which the modified
+ function has a nonpositive function value and nonnegative
+ derivative, and if a lower function value has been
+ obtained but the decrease is not sufficient.
+ */
+ if (stage1 && ftest1 < *f && *f <= fx)
+ {
+ /* Define the modified function and derivative values. */
+ fm = *f - *stp * dgtest;
+ fxm = fx - stx * dgtest;
+ fym = fy - sty * dgtest;
+ dgm = dg - dgtest;
+ dgxm = dgx - dgtest;
+ dgym = dgy - dgtest;
+
+ /*
+ Call update_trial_interval() to update the interval of
+ uncertainty and to compute the new step.
+ */
+ uinfo = update_trial_interval(
+ &stx, &fxm, &dgxm,
+ &sty, &fym, &dgym,
+ stp, &fm, &dgm,
+ stmin, stmax, &brackt);
+
+ /* Reset the function and gradient values for f. */
+ fx = fxm + stx * dgtest;
+ fy = fym + sty * dgtest;
+ dgx = dgxm + dgtest;
+ dgy = dgym + dgtest;
+ }
+ else
+ {
+ /*
+ Call update_trial_interval() to update the interval of
+ uncertainty and to compute the new step.
+ */
+ uinfo = update_trial_interval(
+ &stx, &fx, &dgx,
+ &sty, &fy, &dgy,
+ stp, f, &dg,
+ stmin, stmax, &brackt);
+ }
+
+ /*
+ Force a sufficient decrease in the interval of uncertainty.
+ */
+ if (brackt)
+ {
+ if (0.66 * prev_width <= fabs(sty - stx))
+ {
+ *stp = stx + 0.5 * (sty - stx);
+ }
+ prev_width = width;
+ width = fabs(sty - stx);
+ }
+ }
+
+ return LBFGSERR_LOGICERROR;
+ }
+
+ /**
+ * Default L-BFGS parameters.
+ */
+ static const lbfgs_parameter_t _default_param = {
+ 8,
+ 1e-5,
+ 0,
+ 1e-5,
+ 0,
+ 40,
+ 1e-20,
+ 1e20,
+ 1e-4,
+ 0.9,
+ 1.0e-16,
+ };
+
+ /**
+ * Initialize L-BFGS parameters to the default values.
+ *
+ * Call this function to fill a parameter structure with the default values
+ * and overwrite parameter values if necessary.
+ *
+ * @param param The pointer to the parameter structure.
+ */
+ inline void lbfgs_load_default_parameters(lbfgs_parameter_t *param)
+ {
+ memcpy(param, &_default_param, sizeof(*param));
+ }
+
+ /**
+ * Start a L-BFGS optimization.
+ * A user must implement a function compatible with ::lbfgs_evaluate_t (evaluation
+ * callback) and pass the pointer to the callback function to lbfgs_optimize()
+ * arguments. Similarly, a user can implement a function compatible with
+ * ::lbfgs_stepbound_t to provide an external upper bound for stepsize, and
+ * ::lbfgs_progress_t (progress callback) to obtain the current progress
+ * (e.g., variables, function value, ||G||, etc) and to cancel the iteration
+ * process if necessary. Implementation of the stepbound and the progress callback
+ * is optional: a user can pass NULL if progress notification is not necessary.
+ *
+ * This algorithm terminates an optimization
+ * when:
+ *
+ * ||G|| < g_epsilon \cdot \max(1, ||x||) .
+ *
+ * In this formula, ||.|| denotes the Euclidean norm.
+ *
+ * @param n The number of variables.
+ * @param x The array of variables. A client program can set
+ * default values for the optimization and receive the
+ * optimization result through this array.
+ * @param ptr_fx The pointer to the variable that receives the final
+ * value of the objective function for the variables.
+ * This argument can be set to NULL if the final
+ * value of the objective function is unnecessary.
+ * @param proc_evaluate The callback function to provide function and
+ * gradient evaluations given a current values of
+ * variables. A client program must implement a
+ * callback function compatible with
+ * lbfgs_evaluate_t and pass the pointer to the
+ * callback function.
+ * @param proc_stepbound The callback function to provide values of the
+ * upperbound of the stepsize to search in, provided
+ * with the beginning values of variables before the
+ * linear search, and the current step vector (can
+ * be negative gradient). A client program can implement
+ * this function for more efficient linesearch. If it is
+ * not used, just set it NULL or nullptr.
+ * @param proc_progress The callback function to receive the progress
+ * (the number of iterations, the current value of
+ * the objective function) of the minimization
+ * process. This argument can be set to NULL if
+ * a progress report is unnecessary.
+ * @param instance A user data for the client program. The callback
+ * functions will receive the value of this argument.
+ * @param param The pointer to a structure representing parameters for
+ * L-BFGS optimization. A client program can set this
+ * parameter to NULL to use the default parameters.
+ * Call lbfgs_load_default_parameters() function to
+ * fill a structure with the default values.
+ * @retval int The status code. This function returns zero if the
+ * minimization process terminates without an error. A
+ * non-zero value indicates an error.
+ */
+ inline int lbfgs_optimize(int n,
+ double *x,
+ double *ptr_fx,
+ lbfgs_evaluate_t proc_evaluate,
+ lbfgs_stepbound_t proc_stepbound,
+ lbfgs_progress_t proc_progress,
+ void *instance,
+ lbfgs_parameter_t *_param)
+ {
+ int ret;
+ int i, j, k, ls, end, bound;
+ double step;
+ int loop;
+ double step_min, step_max;
+
+ /* Constant parameters and their default values. */
+ lbfgs_parameter_t param = (_param != NULL) ? (*_param) : _default_param;
+ const int m = param.mem_size;
+
+ double *xp = NULL;
+ double *g = NULL, *gp = NULL;
+ double *d = NULL, *pf = NULL;
+ iteration_data_t *lm = NULL, *it = NULL;
+ double ys, yy;
+ double xnorm, gnorm, beta;
+ double fx = 0.;
+ double rate = 0.;
+
+ /* Construct a callback data. */
+ callback_data_t cd;
+ cd.n = n;
+ cd.instance = instance;
+ cd.proc_evaluate = proc_evaluate;
+ cd.proc_stepbound = proc_stepbound;
+ cd.proc_progress = proc_progress;
+
+ /* Check the input parameters for errors. */
+ if (n <= 0)
+ {
+ return LBFGSERR_INVALID_N;
+ }
+ if (m <= 0)
+ {
+ return LBFGSERR_INVALID_MEMSIZE;
+ }
+ if (param.g_epsilon < 0.)
+ {
+ return LBFGSERR_INVALID_GEPSILON;
+ }
+ if (param.past < 0)
+ {
+ return LBFGSERR_INVALID_TESTPERIOD;
+ }
+ if (param.delta < 0.)
+ {
+ return LBFGSERR_INVALID_DELTA;
+ }
+ if (param.min_step < 0.)
+ {
+ return LBFGSERR_INVALID_MINSTEP;
+ }
+ if (param.max_step < param.min_step)
+ {
+ return LBFGSERR_INVALID_MAXSTEP;
+ }
+ if (param.f_dec_coeff < 0.)
+ {
+ return LBFGSERR_INVALID_FDECCOEFF;
+ }
+ if (param.s_curv_coeff <= param.f_dec_coeff || 1. <= param.s_curv_coeff)
+ {
+ return LBFGSERR_INVALID_SCURVCOEFF;
+ }
+ if (param.xtol < 0.)
+ {
+ return LBFGSERR_INVALID_XTOL;
+ }
+ if (param.max_linesearch <= 0)
+ {
+ return LBFGSERR_INVALID_MAXLINESEARCH;
+ }
+
+ /* Allocate working space. */
+ xp = (double *)vecalloc(n * sizeof(double));
+ g = (double *)vecalloc(n * sizeof(double));
+ gp = (double *)vecalloc(n * sizeof(double));
+ d = (double *)vecalloc(n * sizeof(double));
+
+ /* Allocate limited memory storage. */
+ lm = (iteration_data_t *)vecalloc(m * sizeof(iteration_data_t));
+
+ /* Initialize the limited memory. */
+ for (i = 0; i < m; ++i)
+ {
+ it = &lm[i];
+ it->alpha = 0;
+ it->ys = 0;
+ it->s = (double *)vecalloc(n * sizeof(double));
+ it->y = (double *)vecalloc(n * sizeof(double));
+ }
+
+ /* Allocate an array for storing previous values of the objective function. */
+ if (0 < param.past)
+ {
+ pf = (double *)vecalloc(param.past * sizeof(double));
+ }
+
+ /* Evaluate the function value and its gradient. */
+ fx = cd.proc_evaluate(cd.instance, x, g, cd.n);
+
+ /* Store the initial value of the objective function. */
+ if (pf != NULL)
+ {
+ pf[0] = fx;
+ }
+
+ /*
+ Compute the direction;
+ we assume the initial hessian matrix H_0 as the identity matrix.
+ */
+ vecncpy(d, g, n);
+
+ /*
+ Make sure that the initial variables are not a minimizer.
+ */
+ vec2norm(&xnorm, x, n);
+ vec2norm(&gnorm, g, n);
+
+ if (xnorm < 1.0)
+ xnorm = 1.0;
+ if (gnorm / xnorm <= param.g_epsilon)
+ {
+ ret = LBFGS_ALREADY_MINIMIZED;
+ }
+ else
+ {
+ /* Compute the initial step:
+ step = 1.0 / sqrt(vecdot(d, d, n))
+ */
+ vec2norminv(&step, d, n);
+
+ k = 1;
+ end = 0;
+ loop = 1;
+
+ while (loop == 1)
+ {
+ /* Store the current position and gradient vectors. */
+ veccpy(xp, x, n);
+ veccpy(gp, g, n);
+
+ // If the step bound can be provied dynamically, then apply it.
+ step_min = param.min_step;
+ step_max = param.max_step;
+ if (cd.proc_stepbound)
+ {
+ step_max = cd.proc_stepbound(cd.instance, xp, d, cd.n);
+ step_max = step_max < param.max_step ? step_max : param.max_step;
+ if (step >= step_max)
+ step = step_max / 2.0;
+ }
+
+ /* Search for an optimal step. */
+ ls = line_search_morethuente(n, x, &fx, g, d, &step, xp, gp, &step_min, &step_max, &cd, ¶m);
+
+ if (ls < 0)
+ {
+ /* Revert to the previous point. */
+ veccpy(x, xp, n);
+ veccpy(g, gp, n);
+ ret = ls;
+ loop = 0;
+ continue;
+ }
+
+ /* Compute x and g norms. */
+ vec2norm(&xnorm, x, n);
+ vec2norm(&gnorm, g, n);
+
+ // /* Report the progress. */
+ // if (cd.proc_progress)
+ // {
+ // if ((ret = cd.proc_progress(cd.instance, x, g, fx, xnorm, gnorm, step, cd.n, k, ls)))
+ // {
+ // loop = 0;
+ // continue;
+ // }
+ // }
+
+ /*
+ Convergence test.
+ The criterion is given by the following formula:
+ |g(x)| / \max(1, |x|) < g_epsilon
+ */
+ if (xnorm < 1.0)
+ xnorm = 1.0;
+ if (gnorm / xnorm <= param.g_epsilon)
+ {
+ /* Convergence. */
+ ret = LBFGS_CONVERGENCE;
+ break;
+ }
+
+ /*
+ Test for stopping criterion.
+ The criterion is given by the following formula:
+ |(f(past_x) - f(x))| / f(x) < \delta
+ */
+ if (pf != NULL)
+ {
+ /* We don't test the stopping criterion while k < past. */
+ if (param.past <= k)
+ {
+ /* Compute the relative improvement from the past. */
+ rate = (pf[k % param.past] - fx) / fx;
+
+ /* The stopping criterion. */
+ if (fabs(rate) < param.delta)
+ {
+ ret = LBFGS_STOP;
+ break;
+ }
+ }
+
+ /* Store the current value of the objective function. */
+ pf[k % param.past] = fx;
+ }
+
+ if (param.max_iterations != 0 && param.max_iterations < k + 1)
+ {
+ /* Maximum number of iterations. */
+ ret = LBFGSERR_MAXIMUMITERATION;
+ break;
+ }
+
+ /*
+ Update vectors s and y:
+ s_{k+1} = x_{k+1} - x_{k} = \step * d_{k}.
+ y_{k+1} = g_{k+1} - g_{k}.
+ */
+ it = &lm[end];
+ vecdiff(it->s, x, xp, n);
+ vecdiff(it->y, g, gp, n);
+
+ /*
+ Compute scalars ys and yy:
+ ys = y^t \cdot s = 1 / \rho.
+ yy = y^t \cdot y.
+ Notice that yy is used for scaling the hessian matrix H_0 (Cholesky factor).
+ */
+ vecdot(&ys, it->y, it->s, n);
+ vecdot(&yy, it->y, it->y, n);
+ it->ys = ys;
+
+ /*
+ Recursive formula to compute dir = -(H \cdot g).
+ This is described in page 779 of:
+ Jorge Nocedal.
+ Updating Quasi-Newton Matrices with Limited Storage.
+ Mathematics of Computation, Vol. 35, No. 151,
+ pp. 773--782, 1980.
+ */
+ bound = (m <= k) ? m : k;
+ ++k;
+ end = (end + 1) % m;
+
+ /* Compute the negative of gradients. */
+ vecncpy(d, g, n);
+
+ j = end;
+ for (i = 0; i < bound; ++i)
+ {
+ j = (j + m - 1) % m; /* if (--j == -1) j = m-1; */
+ it = &lm[j];
+ /* \alpha_{j} = \rho_{j} s^{t}_{j} \cdot q_{k+1}. */
+ vecdot(&it->alpha, it->s, d, n);
+ it->alpha /= it->ys;
+ /* q_{i} = q_{i+1} - \alpha_{i} y_{i}. */
+ vecadd(d, it->y, -it->alpha, n);
+ }
+
+ vecscale(d, ys / yy, n);
+
+ for (i = 0; i < bound; ++i)
+ {
+ it = &lm[j];
+ /* \beta_{j} = \rho_{j} y^t_{j} \cdot \gamm_{i}. */
+ vecdot(&beta, it->y, d, n);
+ beta /= it->ys;
+ /* \gamm_{i+1} = \gamm_{i} + (\alpha_{j} - \beta_{j}) s_{j}. */
+ vecadd(d, it->s, it->alpha - beta, n);
+ j = (j + 1) % m; /* if (++j == m) j = 0; */
+ }
+
+ /*
+ Now the search direction d is ready. We try step = 1 first.
+ */
+ step = 1.0;
+ }
+ }
+
+ /* Return the final value of the objective function. */
+ if (ptr_fx != NULL)
+ {
+ *ptr_fx = fx;
+ }
+
+ vecfree(pf);
+
+ /* Free memory blocks used by this function. */
+ if (lm != NULL)
+ {
+ for (i = 0; i < m; ++i)
+ {
+ vecfree(lm[i].s);
+ vecfree(lm[i].y);
+ }
+ vecfree(lm);
+ }
+ vecfree(d);
+ vecfree(gp);
+ vecfree(g);
+ vecfree(xp);
+
+ return ret;
+ }
+
+ /**
+ * Get string description of an lbfgs_optimize() return code.
+ *
+ * @param err A value returned by lbfgs_optimize().
+ */
+ inline const char *lbfgs_strerror(int err)
+ {
+ switch (err)
+ {
+ case LBFGS_CONVERGENCE:
+ return "Success: reached convergence (g_epsilon).";
+
+ case LBFGS_STOP:
+ return "Success: met stopping criteria (past f decrease less than delta).";
+
+ case LBFGS_ALREADY_MINIMIZED:
+ return "The initial variables already minimize the objective function.";
+
+ case LBFGSERR_UNKNOWNERROR:
+ return "Unknown error.";
+
+ case LBFGSERR_LOGICERROR:
+ return "Logic error.";
+
+ case LBFGSERR_CANCELED:
+ return "The minimization process has been canceled.";
+
+ case LBFGSERR_INVALID_N:
+ return "Invalid number of variables specified.";
+
+ case LBFGSERR_INVALID_MEMSIZE:
+ return "Invalid parameter lbfgs_parameter_t::mem_size specified.";
+
+ case LBFGSERR_INVALID_GEPSILON:
+ return "Invalid parameter lbfgs_parameter_t::g_epsilon specified.";
+
+ case LBFGSERR_INVALID_TESTPERIOD:
+ return "Invalid parameter lbfgs_parameter_t::past specified.";
+
+ case LBFGSERR_INVALID_DELTA:
+ return "Invalid parameter lbfgs_parameter_t::delta specified.";
+
+ case LBFGSERR_INVALID_MINSTEP:
+ return "Invalid parameter lbfgs_parameter_t::min_step specified.";
+
+ case LBFGSERR_INVALID_MAXSTEP:
+ return "Invalid parameter lbfgs_parameter_t::max_step specified.";
+
+ case LBFGSERR_INVALID_FDECCOEFF:
+ return "Invalid parameter lbfgs_parameter_t::f_dec_coeff specified.";
+
+ case LBFGSERR_INVALID_SCURVCOEFF:
+ return "Invalid parameter lbfgs_parameter_t::s_curv_coeff specified.";
+
+ case LBFGSERR_INVALID_XTOL:
+ return "Invalid parameter lbfgs_parameter_t::xtol specified.";
+
+ case LBFGSERR_INVALID_MAXLINESEARCH:
+ return "Invalid parameter lbfgs_parameter_t::max_linesearch specified.";
+
+ case LBFGSERR_OUTOFINTERVAL:
+ return "The line-search step went out of the interval of uncertainty.";
+
+ case LBFGSERR_INCORRECT_TMINMAX:
+ return "A logic error occurred; alternatively, the interval of uncertainty"
+ " became too small.";
+
+ case LBFGSERR_ROUNDING_ERROR:
+ return "A rounding error occurred; alternatively, no line-search step"
+ " satisfies the sufficient decrease and curvature conditions.";
+
+ case LBFGSERR_MINIMUMSTEP:
+ return "The line-search step became smaller than lbfgs_parameter_t::min_step.";
+
+ case LBFGSERR_MAXIMUMSTEP:
+ return "The line-search step became larger than lbfgs_parameter_t::max_step.";
+
+ case LBFGSERR_MAXIMUMLINESEARCH:
+ return "The line-search routine reaches the maximum number of evaluations.";
+
+ case LBFGSERR_MAXIMUMITERATION:
+ return "The algorithm routine reaches the maximum number of iterations.";
+
+ case LBFGSERR_WIDTHTOOSMALL:
+ return "Relative width of the interval of uncertainty is at most"
+ " lbfgs_parameter_t::xtol.";
+
+ case LBFGSERR_INVALIDPARAMETERS:
+ return "A logic error (negative line-search step) occurred.";
+
+ case LBFGSERR_INCREASEGRADIENT:
+ return "The current search direction increases the objective function value.";
+
+ default:
+ return "(unknown)";
+ }
+ }
+
+} // namespace lbfgs
+
+#endif
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/uniform_bspline.h b/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/uniform_bspline.h
new file mode 100644
index 0000000..d2507a2
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/bspline_opt/include/bspline_opt/uniform_bspline.h
@@ -0,0 +1,80 @@
+#ifndef _UNIFORM_BSPLINE_H_
+#define _UNIFORM_BSPLINE_H_
+
+#include
+#include
+#include
+
+using namespace std;
+
+namespace ego_planner
+{
+ // An implementation of non-uniform B-spline with different dimensions
+ // It also represents uniform B-spline which is a special case of non-uniform
+ class UniformBspline
+ {
+ private:
+ // control points for B-spline with different dimensions.
+ // Each row represents one single control point
+ // The dimension is determined by column number
+ // e.g. B-spline with N points in 3D space -> Nx3 matrix
+ Eigen::MatrixXd control_points_;
+
+ int p_, n_, m_; // p degree, n+1 control points, m = n+p+1
+ Eigen::VectorXd u_; // knots vector
+ double interval_; // knot span \delta t
+
+ Eigen::MatrixXd getDerivativeControlPoints();
+
+ double limit_vel_, limit_acc_, limit_ratio_, feasibility_tolerance_; // physical limits and time adjustment ratio
+
+ public:
+ UniformBspline() {}
+ UniformBspline(const Eigen::MatrixXd &points, const int &order, const double &interval);
+ ~UniformBspline();
+
+ Eigen::MatrixXd get_control_points(void) { return control_points_; }
+
+ // initialize as an uniform B-spline
+ void setUniformBspline(const Eigen::MatrixXd &points, const int &order, const double &interval);
+
+ // get / set basic bspline info
+
+ void setKnot(const Eigen::VectorXd &knot);
+ Eigen::VectorXd getKnot();
+ Eigen::MatrixXd getControlPoint();
+ double getInterval();
+ bool getTimeSpan(double &um, double &um_p);
+
+ // compute position / derivative
+
+ Eigen::VectorXd evaluateDeBoor(const double &u); // use u \in [up, u_mp]
+ inline Eigen::VectorXd evaluateDeBoorT(const double &t) { return evaluateDeBoor(t + u_(p_)); } // use t \in [0, duration]
+ UniformBspline getDerivative();
+
+ // 3D B-spline interpolation of points in point_set, with boundary vel&acc
+ // constraints
+ // input : (K+2) points with boundary vel/acc; ts
+ // output: (K+6) control_pts
+ static void parameterizeToBspline(const double &ts, const vector &point_set,
+ const vector &start_end_derivative,
+ Eigen::MatrixXd &ctrl_pts);
+
+ /* check feasibility, adjust time */
+
+ void setPhysicalLimits(const double &vel, const double &acc, const double &tolerance);
+ bool checkFeasibility(double &ratio, bool show = false);
+ void lengthenTime(const double &ratio);
+
+ /* for performance evaluation */
+
+ double getTimeSum();
+ double getLength(const double &res = 0.01);
+ double getJerk();
+ void getMeanAndMaxVel(double &mean_v, double &max_v);
+ void getMeanAndMaxAcc(double &mean_a, double &max_a);
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+ };
+} // namespace ego_planner
+#endif
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/bspline_opt/package.xml b/motion_planning/3d/ego_planner/planner/bspline_opt/package.xml
new file mode 100755
index 0000000..02e098b
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/bspline_opt/package.xml
@@ -0,0 +1,77 @@
+
+
+ bspline_opt
+ 0.0.0
+ The bspline_opt package
+
+
+
+
+ iszhouxin
+
+
+
+
+
+ TODO
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ catkin
+ roscpp
+ rospy
+ std_msgs
+ plan_env
+ path_searching
+ traj_utils
+ roscpp
+ rospy
+ std_msgs
+ plan_env
+ path_searching
+ traj_utils
+ roscpp
+ rospy
+ std_msgs
+ plan_env
+ path_searching
+ traj_utils
+
+
+
+
+
+
+
+
diff --git a/motion_planning/3d/ego_planner/planner/bspline_opt/src/bspline_optimizer.cpp b/motion_planning/3d/ego_planner/planner/bspline_opt/src/bspline_optimizer.cpp
new file mode 100644
index 0000000..49cca24
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/bspline_opt/src/bspline_optimizer.cpp
@@ -0,0 +1,1730 @@
+#include "bspline_opt/bspline_optimizer.h"
+#include "bspline_opt/gradient_descent_optimizer.h"
+// using namespace std;
+
+namespace ego_planner
+{
+
+ void BsplineOptimizer::setParam(ros::NodeHandle &nh)
+ {
+ nh.param("optimization/lambda_smooth", lambda1_, -1.0);
+ nh.param("optimization/lambda_collision", lambda2_, -1.0);
+ nh.param("optimization/lambda_feasibility", lambda3_, -1.0);
+ nh.param("optimization/lambda_fitness", lambda4_, -1.0);
+
+ nh.param("optimization/dist0", dist0_, -1.0);
+ nh.param("optimization/swarm_clearance", swarm_clearance_, -1.0);
+ nh.param("optimization/max_vel", max_vel_, -1.0);
+ nh.param("optimization/max_acc", max_acc_, -1.0);
+
+ nh.param("optimization/order", order_, 3);
+ }
+
+ void BsplineOptimizer::setEnvironment(const GridMap::Ptr &map)
+ {
+ this->grid_map_ = map;
+ }
+
+ void BsplineOptimizer::setEnvironment(const GridMap::Ptr &map, const fast_planner::ObjPredictor::Ptr mov_obj)
+ {
+ this->grid_map_ = map;
+ this->moving_objs_ = mov_obj;
+ }
+
+ void BsplineOptimizer::setControlPoints(const Eigen::MatrixXd &points)
+ {
+ cps_.points = points;
+ }
+
+ void BsplineOptimizer::setBsplineInterval(const double &ts) { bspline_interval_ = ts; }
+
+ void BsplineOptimizer::setSwarmTrajs(SwarmTrajData *swarm_trajs_ptr) { swarm_trajs_ = swarm_trajs_ptr; }
+
+ void BsplineOptimizer::setDroneId(const int drone_id) { drone_id_ = drone_id; }
+
+ std::vector BsplineOptimizer::distinctiveTrajs(vector> segments)
+ {
+ if (segments.size() == 0) // will be invoked again later.
+ {
+ std::vector oneSeg;
+ oneSeg.push_back(cps_);
+ return oneSeg;
+ }
+
+ // cout << "A1" << endl;
+
+ constexpr int MAX_TRAJS = 8;
+ constexpr int VARIS = 2;
+ int seg_upbound = std::min((int)segments.size(), static_cast(floor(log(MAX_TRAJS) / log(VARIS))));
+ std::vector control_pts_buf;
+ control_pts_buf.reserve(MAX_TRAJS);
+ const double RESOLUTION = grid_map_->getResolution();
+ const double CTRL_PT_DIST = (cps_.points.col(0) - cps_.points.col(cps_.size - 1)).norm() / (cps_.size - 1);
+
+ // Step 1. Find the opposite vectors and base points for every segment.
+ std::vector> RichInfoSegs;
+ for (int i = 0; i < seg_upbound; i++)
+ {
+ std::pair RichInfoOneSeg;
+ ControlPoints RichInfoOneSeg_temp;
+ cps_.segment(RichInfoOneSeg_temp, segments[i].first, segments[i].second);
+ RichInfoOneSeg.first = RichInfoOneSeg_temp;
+ RichInfoOneSeg.second = RichInfoOneSeg_temp;
+ RichInfoSegs.push_back(RichInfoOneSeg);
+
+ // cout << "RichInfoOneSeg_temp, out" << endl;
+ // cout << "RichInfoSegs[" << i << "].first" << endl;
+ // for ( int k=0; k 0 )
+ // {
+ // cout << "###" << RichInfoOneSeg_temp.points.col(k).transpose() << endl;
+ // for (int k2 = 0; k2 < RichInfoOneSeg_temp.base_point[k].size(); k2++)
+ // {
+ // cout << " " << RichInfoOneSeg_temp.base_point[k][k2].transpose() << " @ " << RichInfoOneSeg_temp.direction[k][k2].transpose() << endl;
+ // }
+ // }
+ }
+
+ // cout << "A2" << endl;
+
+ for (int i = 0; i < seg_upbound; i++)
+ {
+
+ // 1.1 Find the start occupied point id and the last occupied point id
+ if (RichInfoSegs[i].first.size > 1)
+ {
+ int occ_start_id = -1, occ_end_id = -1;
+ Eigen::Vector3d occ_start_pt, occ_end_pt;
+ for (int j = 0; j < RichInfoSegs[i].first.size - 1; j++)
+ {
+ //cout << "A *" << j << "*" << endl;
+ double step_size = RESOLUTION / (RichInfoSegs[i].first.points.col(j) - RichInfoSegs[i].first.points.col(j + 1)).norm() / 2;
+ for (double a = 1; a > 0; a -= step_size)
+ {
+ Eigen::Vector3d pt(a * RichInfoSegs[i].first.points.col(j) + (1 - a) * RichInfoSegs[i].first.points.col(j + 1));
+ //cout << " " << grid_map_->getInflateOccupancy(pt) << " pt=" << pt.transpose() << endl;
+ if (grid_map_->getInflateOccupancy(pt))
+ {
+ occ_start_id = j;
+ occ_start_pt = pt;
+ goto exit_multi_loop1;
+ }
+ }
+ }
+ exit_multi_loop1:;
+ for (int j = RichInfoSegs[i].first.size - 1; j >= 1; j--)
+ {
+ //cout << "j=" << j << endl;
+ //cout << "B *" << j << "*" << endl;
+ ;
+ double step_size = RESOLUTION / (RichInfoSegs[i].first.points.col(j) - RichInfoSegs[i].first.points.col(j - 1)).norm();
+ for (double a = 1; a > 0; a -= step_size)
+ {
+ Eigen::Vector3d pt(a * RichInfoSegs[i].first.points.col(j) + (1 - a) * RichInfoSegs[i].first.points.col(j - 1));
+ //cout << " " << grid_map_->getInflateOccupancy(pt) << " pt=" << pt.transpose() << endl;
+ ;
+ if (grid_map_->getInflateOccupancy(pt))
+ {
+ occ_end_id = j;
+ occ_end_pt = pt;
+ goto exit_multi_loop2;
+ }
+ }
+ }
+ exit_multi_loop2:;
+
+ // double check
+ if (occ_start_id == -1 || occ_end_id == -1)
+ {
+ // It means that the first or the last control points of one segment are in obstacles, which is not allowed.
+ // ROS_WARN("What? occ_start_id=%d, occ_end_id=%d", occ_start_id, occ_end_id);
+
+ segments.erase(segments.begin() + i);
+ RichInfoSegs.erase(RichInfoSegs.begin() + i);
+ seg_upbound--;
+ i--;
+
+ continue;
+
+ // cout << "RichInfoSegs[" << i << "].first" << endl;
+ // for (int k = 0; k < RichInfoSegs[i].first.size; k++)
+ // {
+ // if (RichInfoSegs[i].first.base_point.size() > 0)
+ // {
+ // cout << "###" << RichInfoSegs[i].first.points.col(k).transpose() << endl;
+ // for (int k2 = 0; k2 < RichInfoSegs[i].first.base_point[k].size(); k2++)
+ // {
+ // cout << " " << RichInfoSegs[i].first.base_point[k][k2].transpose() << " @ " << RichInfoSegs[i].first.direction[k][k2].transpose() << endl;
+ // }
+ // }
+ // }
+ }
+
+ // 1.2 Reverse the vector and find new base points from occ_start_id to occ_end_id.
+ for (int j = occ_start_id; j <= occ_end_id; j++)
+ {
+ Eigen::Vector3d base_pt_reverse, base_vec_reverse;
+ if (RichInfoSegs[i].first.base_point[j].size() != 1)
+ {
+ cout << "RichInfoSegs[" << i << "].first.base_point[" << j << "].size()=" << RichInfoSegs[i].first.base_point[j].size() << endl;
+ ROS_ERROR("Wrong number of base_points!!! Should not be happen!.");
+
+ cout << setprecision(5);
+ cout << "cps_" << endl;
+ cout << " clearance=" << cps_.clearance << " cps.size=" << cps_.size << endl;
+ for (int temp_i = 0; temp_i < cps_.size; temp_i++)
+ {
+ if (cps_.base_point[temp_i].size() > 1 && cps_.base_point[temp_i].size() < 1000)
+ {
+ ROS_ERROR("Should not happen!!!");
+ cout << "######" << cps_.points.col(temp_i).transpose() << endl;
+ for (size_t temp_j = 0; temp_j < cps_.base_point[temp_i].size(); temp_j++)
+ cout << " " << cps_.base_point[temp_i][temp_j].transpose() << " @ " << cps_.direction[temp_i][temp_j].transpose() << endl;
+ }
+ }
+
+ std::vector blank;
+ return blank;
+ }
+
+ base_vec_reverse = -RichInfoSegs[i].first.direction[j][0];
+
+ // The start and the end case must get taken special care of.
+ if (j == occ_start_id)
+ {
+ base_pt_reverse = occ_start_pt;
+ }
+ else if (j == occ_end_id)
+ {
+ base_pt_reverse = occ_end_pt;
+ }
+ else
+ {
+ base_pt_reverse = RichInfoSegs[i].first.points.col(j) + base_vec_reverse * (RichInfoSegs[i].first.base_point[j][0] - RichInfoSegs[i].first.points.col(j)).norm();
+ }
+
+ if (grid_map_->getInflateOccupancy(base_pt_reverse)) // Search outward.
+ {
+ double l_upbound = 5 * CTRL_PT_DIST; // "5" is the threshold.
+ double l = RESOLUTION;
+ for (; l <= l_upbound; l += RESOLUTION)
+ {
+ Eigen::Vector3d base_pt_temp = base_pt_reverse + l * base_vec_reverse;
+ //cout << base_pt_temp.transpose() << endl;
+ if (!grid_map_->getInflateOccupancy(base_pt_temp))
+ {
+ RichInfoSegs[i].second.base_point[j][0] = base_pt_temp;
+ RichInfoSegs[i].second.direction[j][0] = base_vec_reverse;
+ break;
+ }
+ }
+ if (l > l_upbound)
+ {
+ ROS_WARN("Can't find the new base points at the opposite within the threshold. i=%d, j=%d", i, j);
+
+ segments.erase(segments.begin() + i);
+ RichInfoSegs.erase(RichInfoSegs.begin() + i);
+ seg_upbound--;
+ i--;
+
+ goto exit_multi_loop3; // break "for (int j = 0; j < RichInfoSegs[i].first.size; j++)"
+ }
+ }
+ else if ((base_pt_reverse - RichInfoSegs[i].first.points.col(j)).norm() >= RESOLUTION) // Unnecessary to search.
+ {
+ RichInfoSegs[i].second.base_point[j][0] = base_pt_reverse;
+ RichInfoSegs[i].second.direction[j][0] = base_vec_reverse;
+ }
+ else
+ {
+ ROS_WARN("base_point and control point are too close!");
+ cout << "base_point=" << RichInfoSegs[i].first.base_point[j][0].transpose() << " control point=" << RichInfoSegs[i].first.points.col(j).transpose() << endl;
+
+ segments.erase(segments.begin() + i);
+ RichInfoSegs.erase(RichInfoSegs.begin() + i);
+ seg_upbound--;
+ i--;
+
+ goto exit_multi_loop3; // break "for (int j = 0; j < RichInfoSegs[i].first.size; j++)"
+ }
+ }
+
+ // 1.3 Assign the base points to control points within [0, occ_start_id) and (occ_end_id, RichInfoSegs[i].first.size()-1].
+ if (RichInfoSegs[i].second.size)
+ {
+ for (int j = occ_start_id - 1; j >= 0; j--)
+ {
+ RichInfoSegs[i].second.base_point[j][0] = RichInfoSegs[i].second.base_point[occ_start_id][0];
+ RichInfoSegs[i].second.direction[j][0] = RichInfoSegs[i].second.direction[occ_start_id][0];
+ }
+ for (int j = occ_end_id + 1; j < RichInfoSegs[i].second.size; j++)
+ {
+ RichInfoSegs[i].second.base_point[j][0] = RichInfoSegs[i].second.base_point[occ_end_id][0];
+ RichInfoSegs[i].second.direction[j][0] = RichInfoSegs[i].second.direction[occ_end_id][0];
+ }
+ }
+
+ exit_multi_loop3:;
+ }
+ else
+ {
+ Eigen::Vector3d base_vec_reverse = -RichInfoSegs[i].first.direction[0][0];
+ Eigen::Vector3d base_pt_reverse = RichInfoSegs[i].first.points.col(0) + base_vec_reverse * (RichInfoSegs[i].first.base_point[0][0] - RichInfoSegs[i].first.points.col(0)).norm();
+
+ if (grid_map_->getInflateOccupancy(base_pt_reverse)) // Search outward.
+ {
+ double l_upbound = 5 * CTRL_PT_DIST; // "5" is the threshold.
+ double l = RESOLUTION;
+ for (; l <= l_upbound; l += RESOLUTION)
+ {
+ Eigen::Vector3d base_pt_temp = base_pt_reverse + l * base_vec_reverse;
+ //cout << base_pt_temp.transpose() << endl;
+ if (!grid_map_->getInflateOccupancy(base_pt_temp))
+ {
+ RichInfoSegs[i].second.base_point[0][0] = base_pt_temp;
+ RichInfoSegs[i].second.direction[0][0] = base_vec_reverse;
+ break;
+ }
+ }
+ if (l > l_upbound)
+ {
+ ROS_WARN("Can't find the new base points at the opposite within the threshold, 2. i=%d", i);
+
+ segments.erase(segments.begin() + i);
+ RichInfoSegs.erase(RichInfoSegs.begin() + i);
+ seg_upbound--;
+ i--;
+ }
+ }
+ else if ((base_pt_reverse - RichInfoSegs[i].first.points.col(0)).norm() >= RESOLUTION) // Unnecessary to search.
+ {
+ RichInfoSegs[i].second.base_point[0][0] = base_pt_reverse;
+ RichInfoSegs[i].second.direction[0][0] = base_vec_reverse;
+ }
+ else
+ {
+ ROS_WARN("base_point and control point are too close!, 2");
+ cout << "base_point=" << RichInfoSegs[i].first.base_point[0][0].transpose() << " control point=" << RichInfoSegs[i].first.points.col(0).transpose() << endl;
+
+ segments.erase(segments.begin() + i);
+ RichInfoSegs.erase(RichInfoSegs.begin() + i);
+ seg_upbound--;
+ i--;
+ }
+ }
+
+ // cout << "A3" << endl;
+ }
+
+ // Step 2. Assemble each segment to make up the new control point sequence.
+ if (seg_upbound == 0) // After the erase operation above, segment legth will decrease to 0 again.
+ {
+ std::vector oneSeg;
+ oneSeg.push_back(cps_);
+ return oneSeg;
+ }
+
+ // cout << "A4" << endl;
+
+ std::vector selection(seg_upbound);
+ std::fill(selection.begin(), selection.end(), 0);
+ selection[0] = -1; // init
+ int max_traj_nums = static_cast(pow(VARIS, seg_upbound));
+ for (int i = 0; i < max_traj_nums; i++)
+ {
+ // 2.1 Calculate the selection table.
+ int digit_id = 0;
+ selection[digit_id]++;
+ while (digit_id < seg_upbound && selection[digit_id] >= VARIS)
+ {
+ selection[digit_id] = 0;
+ digit_id++;
+ if (digit_id >= seg_upbound)
+ {
+ ROS_ERROR("Should not happen!!! digit_id=%d, seg_upbound=%d", digit_id, seg_upbound);
+ }
+ selection[digit_id]++;
+ }
+
+ // 2.2 Assign params according to the selection table.
+ ControlPoints cpsOneSample;
+ cpsOneSample.resize(cps_.size);
+ cpsOneSample.clearance = cps_.clearance;
+ int cp_id = 0, seg_id = 0, cp_of_seg_id = 0;
+ while (/*seg_id < RichInfoSegs.size() ||*/ cp_id < cps_.size)
+ {
+ //cout << "A ";
+ // if ( seg_id >= RichInfoSegs.size() )
+ // {
+ // cout << "seg_id=" << seg_id << " RichInfoSegs.size()=" << RichInfoSegs.size() << endl;
+ // }
+ // if ( cp_id >= cps_.base_point.size() )
+ // {
+ // cout << "cp_id=" << cp_id << " cps_.base_point.size()=" << cps_.base_point.size() << endl;
+ // }
+ // if ( cp_of_seg_id >= RichInfoSegs[seg_id].first.base_point.size() )
+ // {
+ // cout << "cp_of_seg_id=" << cp_of_seg_id << " RichInfoSegs[seg_id].first.base_point.size()=" << RichInfoSegs[seg_id].first.base_point.size() << endl;
+ // }
+
+ if (seg_id >= seg_upbound || cp_id < segments[seg_id].first || cp_id > segments[seg_id].second)
+ {
+ cpsOneSample.points.col(cp_id) = cps_.points.col(cp_id);
+ cpsOneSample.base_point[cp_id] = cps_.base_point[cp_id];
+ cpsOneSample.direction[cp_id] = cps_.direction[cp_id];
+ }
+ else if (cp_id >= segments[seg_id].first && cp_id <= segments[seg_id].second)
+ {
+ if (!selection[seg_id]) // zx-todo
+ {
+ cpsOneSample.points.col(cp_id) = RichInfoSegs[seg_id].first.points.col(cp_of_seg_id);
+ cpsOneSample.base_point[cp_id] = RichInfoSegs[seg_id].first.base_point[cp_of_seg_id];
+ cpsOneSample.direction[cp_id] = RichInfoSegs[seg_id].first.direction[cp_of_seg_id];
+ cp_of_seg_id++;
+ }
+ else
+ {
+ if (RichInfoSegs[seg_id].second.size)
+ {
+ cpsOneSample.points.col(cp_id) = RichInfoSegs[seg_id].second.points.col(cp_of_seg_id);
+ cpsOneSample.base_point[cp_id] = RichInfoSegs[seg_id].second.base_point[cp_of_seg_id];
+ cpsOneSample.direction[cp_id] = RichInfoSegs[seg_id].second.direction[cp_of_seg_id];
+ cp_of_seg_id++;
+ }
+ else
+ {
+ // Abandon this trajectory.
+ goto abandon_this_trajectory;
+ }
+ }
+
+ if (cp_id == segments[seg_id].second)
+ {
+ cp_of_seg_id = 0;
+ seg_id++;
+ }
+ }
+ else
+ {
+ ROS_ERROR("Shold not happen!!!!, cp_id=%d, seg_id=%d, segments.front().first=%d, segments.back().second=%d, segments[seg_id].first=%d, segments[seg_id].second=%d",
+ cp_id, seg_id, segments.front().first, segments.back().second, segments[seg_id].first, segments[seg_id].second);
+ }
+
+ cp_id++;
+ }
+
+ control_pts_buf.push_back(cpsOneSample);
+
+ abandon_this_trajectory:;
+ }
+
+ // cout << "A5" << endl;
+
+ return control_pts_buf;
+ } // namespace ego_planner
+
+ /* This function is very similar to check_collision_and_rebound().
+ * It was written separately, just because I did it once and it has been running stably since March 2020.
+ * But I will merge then someday.*/
+ std::vector> BsplineOptimizer::initControlPoints(Eigen::MatrixXd &init_points, bool flag_first_init /*= true*/)
+ {
+
+ if (flag_first_init)
+ {
+ cps_.clearance = dist0_;
+ cps_.resize(init_points.cols());
+ cps_.points = init_points;
+ }
+
+ /*** Segment the initial trajectory according to obstacles ***/
+ constexpr int ENOUGH_INTERVAL = 2;
+ double step_size = grid_map_->getResolution() / ((init_points.col(0) - init_points.rightCols(1)).norm() / (init_points.cols() - 1)) / 1.5;
+ int in_id, out_id;
+ vector> segment_ids;
+ int same_occ_state_times = ENOUGH_INTERVAL + 1;
+ bool occ, last_occ = false;
+ bool flag_got_start = false, flag_got_end = false, flag_got_end_maybe = false;
+ int i_end = (int)init_points.cols() - order_ - ((int)init_points.cols() - 2 * order_) / 3; // only check closed 2/3 points.
+ for (int i = order_; i <= i_end; ++i)
+ {
+ //cout << " *" << i-1 << "*" ;
+ for (double a = 1.0; a > 0.0; a -= step_size)
+ {
+ occ = grid_map_->getInflateOccupancy(a * init_points.col(i - 1) + (1 - a) * init_points.col(i));
+ //cout << " " << occ;
+ // cout << setprecision(5);
+ // cout << (a * init_points.col(i-1) + (1-a) * init_points.col(i)).transpose() << " occ1=" << occ << endl;
+
+ if (occ && !last_occ)
+ {
+ if (same_occ_state_times > ENOUGH_INTERVAL || i == order_)
+ {
+ in_id = i - 1;
+ flag_got_start = true;
+ }
+ same_occ_state_times = 0;
+ flag_got_end_maybe = false; // terminate in advance
+ }
+ else if (!occ && last_occ)
+ {
+ out_id = i;
+ flag_got_end_maybe = true;
+ same_occ_state_times = 0;
+ }
+ else
+ {
+ ++same_occ_state_times;
+ }
+
+ if (flag_got_end_maybe && (same_occ_state_times > ENOUGH_INTERVAL || (i == (int)init_points.cols() - order_)))
+ {
+ flag_got_end_maybe = false;
+ flag_got_end = true;
+ }
+
+ last_occ = occ;
+
+ if (flag_got_start && flag_got_end)
+ {
+ flag_got_start = false;
+ flag_got_end = false;
+ segment_ids.push_back(std::pair(in_id, out_id));
+ }
+ }
+ }
+ // cout << endl;
+
+ // for (size_t i = 0; i < segment_ids.size(); i++)
+ // {
+ // cout << "segment_ids=" << segment_ids[i].first << " ~ " << segment_ids[i].second << endl;
+ // }
+
+ // return in advance
+ if (segment_ids.size() == 0)
+ {
+ vector> blank_ret;
+ return blank_ret;
+ }
+
+ /*** a star search ***/
+ vector> a_star_pathes;
+ for (size_t i = 0; i < segment_ids.size(); ++i)
+ {
+ //cout << "in=" << in.transpose() << " out=" << out.transpose() << endl;
+ Eigen::Vector3d in(init_points.col(segment_ids[i].first)), out(init_points.col(segment_ids[i].second));
+ if (a_star_->AstarSearch(/*(in-out).norm()/10+0.05*/ 0.1, in, out))
+ {
+ a_star_pathes.push_back(a_star_->getPath());
+ }
+ else
+ {
+ ROS_ERROR("a star error, force return!");
+ vector> blank_ret;
+ return blank_ret;
+ }
+ }
+
+ /*** calculate bounds ***/
+ int id_low_bound, id_up_bound;
+ vector> bounds(segment_ids.size());
+ for (size_t i = 0; i < segment_ids.size(); i++)
+ {
+
+ if (i == 0) // first segment
+ {
+ id_low_bound = order_;
+ if (segment_ids.size() > 1)
+ {
+ id_up_bound = (int)(((segment_ids[0].second + segment_ids[1].first) - 1.0f) / 2); // id_up_bound : -1.0f fix()
+ }
+ else
+ {
+ id_up_bound = init_points.cols() - order_ - 1;
+ }
+ }
+ else if (i == segment_ids.size() - 1) // last segment, i != 0 here
+ {
+ id_low_bound = (int)(((segment_ids[i].first + segment_ids[i - 1].second) + 1.0f) / 2); // id_low_bound : +1.0f ceil()
+ id_up_bound = init_points.cols() - order_ - 1;
+ }
+ else
+ {
+ id_low_bound = (int)(((segment_ids[i].first + segment_ids[i - 1].second) + 1.0f) / 2); // id_low_bound : +1.0f ceil()
+ id_up_bound = (int)(((segment_ids[i].second + segment_ids[i + 1].first) - 1.0f) / 2); // id_up_bound : -1.0f fix()
+ }
+
+ bounds[i] = std::pair(id_low_bound, id_up_bound);
+ }
+
+ // cout << "+++++++++" << endl;
+ // for ( int j=0; j> adjusted_segment_ids(segment_ids.size());
+ constexpr double MINIMUM_PERCENT = 0.0; // Each segment is guaranteed to have sufficient points to generate sufficient force
+ int minimum_points = round(init_points.cols() * MINIMUM_PERCENT), num_points;
+ for (size_t i = 0; i < segment_ids.size(); i++)
+ {
+ /*** Adjust segment length ***/
+ num_points = segment_ids[i].second - segment_ids[i].first + 1;
+ //cout << "i = " << i << " first = " << segment_ids[i].first << " second = " << segment_ids[i].second << endl;
+ if (num_points < minimum_points)
+ {
+ double add_points_each_side = (int)(((minimum_points - num_points) + 1.0f) / 2);
+
+ adjusted_segment_ids[i].first = segment_ids[i].first - add_points_each_side >= bounds[i].first ? segment_ids[i].first - add_points_each_side : bounds[i].first;
+
+ adjusted_segment_ids[i].second = segment_ids[i].second + add_points_each_side <= bounds[i].second ? segment_ids[i].second + add_points_each_side : bounds[i].second;
+ }
+ else
+ {
+ adjusted_segment_ids[i].first = segment_ids[i].first;
+ adjusted_segment_ids[i].second = segment_ids[i].second;
+ }
+
+ //cout << "final:" << "i = " << i << " first = " << adjusted_segment_ids[i].first << " second = " << adjusted_segment_ids[i].second << endl;
+ }
+ for (size_t i = 1; i < adjusted_segment_ids.size(); i++) // Avoid overlap
+ {
+ if (adjusted_segment_ids[i - 1].second >= adjusted_segment_ids[i].first)
+ {
+ double middle = (double)(adjusted_segment_ids[i - 1].second + adjusted_segment_ids[i].first) / 2.0;
+ adjusted_segment_ids[i - 1].second = static_cast(middle - 0.1);
+ adjusted_segment_ids[i].first = static_cast(middle + 1.1);
+ }
+ }
+
+ // Used for return
+ vector> final_segment_ids;
+
+ /*** Assign data to each segment ***/
+ for (size_t i = 0; i < segment_ids.size(); i++)
+ {
+ // step 1
+ for (int j = adjusted_segment_ids[i].first; j <= adjusted_segment_ids[i].second; ++j)
+ cps_.flag_temp[j] = false;
+
+ // step 2
+ int got_intersection_id = -1;
+ for (int j = segment_ids[i].first + 1; j < segment_ids[i].second; ++j)
+ {
+ Eigen::Vector3d ctrl_pts_law(init_points.col(j + 1) - init_points.col(j - 1)), intersection_point;
+ int Astar_id = a_star_pathes[i].size() / 2, last_Astar_id; // Let "Astar_id = id_of_the_most_far_away_Astar_point" will be better, but it needs more computation
+ double val = (a_star_pathes[i][Astar_id] - init_points.col(j)).dot(ctrl_pts_law), last_val = val;
+ while (Astar_id >= 0 && Astar_id < (int)a_star_pathes[i].size())
+ {
+ last_Astar_id = Astar_id;
+
+ if (val >= 0)
+ --Astar_id;
+ else
+ ++Astar_id;
+
+ val = (a_star_pathes[i][Astar_id] - init_points.col(j)).dot(ctrl_pts_law);
+
+ if (val * last_val <= 0 && (abs(val) > 0 || abs(last_val) > 0)) // val = last_val = 0.0 is not allowed
+ {
+ intersection_point =
+ a_star_pathes[i][Astar_id] +
+ ((a_star_pathes[i][Astar_id] - a_star_pathes[i][last_Astar_id]) *
+ (ctrl_pts_law.dot(init_points.col(j) - a_star_pathes[i][Astar_id]) / ctrl_pts_law.dot(a_star_pathes[i][Astar_id] - a_star_pathes[i][last_Astar_id])) // = t
+ );
+
+ //cout << "i=" << i << " j=" << j << " Astar_id=" << Astar_id << " last_Astar_id=" << last_Astar_id << " intersection_point = " << intersection_point.transpose() << endl;
+
+ got_intersection_id = j;
+ break;
+ }
+ }
+
+ if (got_intersection_id >= 0)
+ {
+ double length = (intersection_point - init_points.col(j)).norm();
+ if (length > 1e-5)
+ {
+ cps_.flag_temp[j] = true;
+ for (double a = length; a >= 0.0; a -= grid_map_->getResolution())
+ {
+ occ = grid_map_->getInflateOccupancy((a / length) * intersection_point + (1 - a / length) * init_points.col(j));
+
+ if (occ || a < grid_map_->getResolution())
+ {
+ if (occ)
+ a += grid_map_->getResolution();
+ cps_.base_point[j].push_back((a / length) * intersection_point + (1 - a / length) * init_points.col(j));
+ cps_.direction[j].push_back((intersection_point - init_points.col(j)).normalized());
+ // cout << "A " << j << endl;
+ break;
+ }
+ }
+ }
+ else
+ {
+ got_intersection_id = -1;
+ }
+ }
+ }
+
+ /* Corner case: the segment length is too short. Here the control points may outside the A* path, leading to opposite gradient direction. So I have to take special care of it */
+ if (segment_ids[i].second - segment_ids[i].first == 1)
+ {
+ Eigen::Vector3d ctrl_pts_law(init_points.col(segment_ids[i].second) - init_points.col(segment_ids[i].first)), intersection_point;
+ Eigen::Vector3d middle_point = (init_points.col(segment_ids[i].second) + init_points.col(segment_ids[i].first)) / 2;
+ int Astar_id = a_star_pathes[i].size() / 2, last_Astar_id; // Let "Astar_id = id_of_the_most_far_away_Astar_point" will be better, but it needs more computation
+ double val = (a_star_pathes[i][Astar_id] - middle_point).dot(ctrl_pts_law), last_val = val;
+ while (Astar_id >= 0 && Astar_id < (int)a_star_pathes[i].size())
+ {
+ last_Astar_id = Astar_id;
+
+ if (val >= 0)
+ --Astar_id;
+ else
+ ++Astar_id;
+
+ val = (a_star_pathes[i][Astar_id] - middle_point).dot(ctrl_pts_law);
+
+ if (val * last_val <= 0 && (abs(val) > 0 || abs(last_val) > 0)) // val = last_val = 0.0 is not allowed
+ {
+ intersection_point =
+ a_star_pathes[i][Astar_id] +
+ ((a_star_pathes[i][Astar_id] - a_star_pathes[i][last_Astar_id]) *
+ (ctrl_pts_law.dot(middle_point - a_star_pathes[i][Astar_id]) / ctrl_pts_law.dot(a_star_pathes[i][Astar_id] - a_star_pathes[i][last_Astar_id])) // = t
+ );
+
+ if ((intersection_point - middle_point).norm() > 0.01) // 1cm.
+ {
+ cps_.flag_temp[segment_ids[i].first] = true;
+ cps_.base_point[segment_ids[i].first].push_back(init_points.col(segment_ids[i].first));
+ cps_.direction[segment_ids[i].first].push_back((intersection_point - middle_point).normalized());
+ // cout << "AA " << segment_ids[i].first << endl;
+
+ got_intersection_id = segment_ids[i].first;
+ }
+ break;
+ }
+ }
+ }
+
+ //step 3
+ if (got_intersection_id >= 0)
+ {
+ for (int j = got_intersection_id + 1; j <= adjusted_segment_ids[i].second; ++j)
+ if (!cps_.flag_temp[j])
+ {
+ cps_.base_point[j].push_back(cps_.base_point[j - 1].back());
+ cps_.direction[j].push_back(cps_.direction[j - 1].back());
+ // cout << "AAA " << j << endl;
+ }
+
+ for (int j = got_intersection_id - 1; j >= adjusted_segment_ids[i].first; --j)
+ if (!cps_.flag_temp[j])
+ {
+ cps_.base_point[j].push_back(cps_.base_point[j + 1].back());
+ cps_.direction[j].push_back(cps_.direction[j + 1].back());
+ // cout << "AAAA " << j << endl;
+ }
+
+ final_segment_ids.push_back(adjusted_segment_ids[i]);
+ }
+ else
+ {
+ // Just ignore, it does not matter ^_^.
+ // ROS_ERROR("Failed to generate direction! segment_id=%d", i);
+ }
+ }
+
+ return final_segment_ids;
+ }
+
+ int BsplineOptimizer::earlyExit(void *func_data, const double *x, const double *g, const double fx, const double xnorm, const double gnorm, const double step, int n, int k, int ls)
+ {
+ BsplineOptimizer *opt = reinterpret_cast(func_data);
+ // cout << "k=" << k << endl;
+ // cout << "opt->flag_continue_to_optimize_=" << opt->flag_continue_to_optimize_ << endl;
+ return (opt->force_stop_type_ == STOP_FOR_ERROR || opt->force_stop_type_ == STOP_FOR_REBOUND);
+ }
+
+ double BsplineOptimizer::costFunctionRebound(void *func_data, const double *x, double *grad, const int n)
+ {
+ BsplineOptimizer *opt = reinterpret_cast(func_data);
+
+ double cost;
+ opt->combineCostRebound(x, grad, cost, n);
+
+ opt->iter_num_ += 1;
+ return cost;
+ }
+
+ double BsplineOptimizer::costFunctionRefine(void *func_data, const double *x, double *grad, const int n)
+ {
+ BsplineOptimizer *opt = reinterpret_cast(func_data);
+
+ double cost;
+ opt->combineCostRefine(x, grad, cost, n);
+
+ opt->iter_num_ += 1;
+ return cost;
+ }
+
+ void BsplineOptimizer::calcSwarmCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient)
+ {
+ cost = 0.0;
+ int end_idx = q.cols() - order_ - (double)(q.cols() - 2*order_)*1.0/3.0; // Only check the first 2/3 points
+ const double CLEARANCE = swarm_clearance_ * 2;
+ double t_now = ros::Time::now().toSec();
+ constexpr double a = 2.0, b = 1.0, inv_a2 = 1 / a / a, inv_b2 = 1 / b / b;
+
+ for (int i = order_; i < end_idx; i++)
+ {
+ double glb_time = t_now + ((double)(order_ - 1) / 2 + (i - order_ + 1)) * bspline_interval_;
+
+ for (size_t id = 0; id < swarm_trajs_->size(); id++)
+ {
+ if ( (swarm_trajs_->at(id).drone_id != (int)id) || swarm_trajs_->at(id).drone_id == drone_id_ )
+ {
+ continue;
+ }
+
+ double traj_i_satrt_time = swarm_trajs_->at(id).start_time_.toSec();
+ if (glb_time < traj_i_satrt_time + swarm_trajs_->at(id).duration_ - 0.1)
+ {
+ /* def cost=(c-sqrt([Q-O]'D[Q-O]))^2, D=[1/b^2,0,0;0,1/b^2,0;0,0,1/a^2] */
+ Eigen::Vector3d swarm_prid = swarm_trajs_->at(id).position_traj_.evaluateDeBoorT(glb_time - traj_i_satrt_time);
+ Eigen::Vector3d dist_vec = cps_.points.col(i) - swarm_prid;
+ double ellip_dist = sqrt(dist_vec(2) * dist_vec(2) * inv_a2 + (dist_vec(0) * dist_vec(0) + dist_vec(1) * dist_vec(1)) * inv_b2);
+ double dist_err = CLEARANCE - ellip_dist;
+
+ Eigen::Vector3d dist_grad = cps_.points.col(i) - swarm_prid;
+ Eigen::Vector3d Coeff;
+ Coeff(0) = -2 * (CLEARANCE / ellip_dist - 1) * inv_b2;
+ Coeff(1) = Coeff(0);
+ Coeff(2) = -2 * (CLEARANCE / ellip_dist - 1) * inv_a2;
+
+ if (dist_err < 0)
+ {
+ /* do nothing */
+ }
+ else
+ {
+ cost += pow(dist_err, 2);
+ gradient.col(i) += (Coeff.array() * dist_grad.array()).matrix();
+ }
+
+ if (min_ellip_dist_ > dist_err)
+ {
+ min_ellip_dist_ = dist_err;
+ }
+ }
+ }
+ }
+ }
+
+ void BsplineOptimizer::calcMovingObjCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient)
+ {
+ cost = 0.0;
+ int end_idx = q.cols() - order_;
+ constexpr double CLEARANCE = 1.5;
+ double t_now = ros::Time::now().toSec();
+
+ for (int i = order_; i < end_idx; i++)
+ {
+ double time = ((double)(order_ - 1) / 2 + (i - order_ + 1)) * bspline_interval_;
+
+ for (int id = 0; id < moving_objs_->getObjNums(); id++)
+ {
+ Eigen::Vector3d obj_prid = moving_objs_->evaluateConstVel(id, t_now + time);
+ double dist = (cps_.points.col(i) - obj_prid).norm();
+ //cout /*<< "cps_.points.col(i)=" << cps_.points.col(i).transpose()*/ << " moving_objs_=" << obj_prid.transpose() << " dist=" << dist << endl;
+ double dist_err = CLEARANCE - dist;
+ Eigen::Vector3d dist_grad = (cps_.points.col(i) - obj_prid).normalized();
+
+ if (dist_err < 0)
+ {
+ /* do nothing */
+ }
+ else
+ {
+ cost += pow(dist_err, 2);
+ gradient.col(i) += -2.0 * dist_err * dist_grad;
+ }
+ }
+ // cout << "time=" << time << " i=" << i << " order_=" << order_ << " end_idx=" << end_idx << endl;
+ // cout << "--" << endl;
+ }
+ // cout << "---------------" << endl;
+ }
+
+ void BsplineOptimizer::calcDistanceCostRebound(const Eigen::MatrixXd &q, double &cost,
+ Eigen::MatrixXd &gradient, int iter_num, double smoothness_cost)
+ {
+ cost = 0.0;
+ int end_idx = q.cols() - order_;
+ double demarcation = cps_.clearance;
+ double a = 3 * demarcation, b = -3 * pow(demarcation, 2), c = pow(demarcation, 3);
+
+ force_stop_type_ = DONT_STOP;
+ if (iter_num > 3 && smoothness_cost / (cps_.size - 2 * order_) < 0.1) // 0.1 is an experimental value that indicates the trajectory is smooth enough.
+ {
+ check_collision_and_rebound();
+ }
+
+ /*** calculate distance cost and gradient ***/
+ for (auto i = order_; i < end_idx; ++i)
+ {
+ for (size_t j = 0; j < cps_.direction[i].size(); ++j)
+ {
+ double dist = (cps_.points.col(i) - cps_.base_point[i][j]).dot(cps_.direction[i][j]);
+ double dist_err = cps_.clearance - dist;
+ Eigen::Vector3d dist_grad = cps_.direction[i][j];
+
+ if (dist_err < 0)
+ {
+ /* do nothing */
+ }
+ else if (dist_err < demarcation)
+ {
+ cost += pow(dist_err, 3);
+ gradient.col(i) += -3.0 * dist_err * dist_err * dist_grad;
+ }
+ else
+ {
+ cost += a * dist_err * dist_err + b * dist_err + c;
+ gradient.col(i) += -(2.0 * a * dist_err + b) * dist_grad;
+ }
+ }
+ }
+ }
+
+ void BsplineOptimizer::calcFitnessCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient)
+ {
+
+ cost = 0.0;
+
+ int end_idx = q.cols() - order_;
+
+ // def: f = |x*v|^2/a^2 + |x×v|^2/b^2
+ double a2 = 25, b2 = 1;
+ for (auto i = order_ - 1; i < end_idx + 1; ++i)
+ {
+ Eigen::Vector3d x = (q.col(i - 1) + 4 * q.col(i) + q.col(i + 1)) / 6.0 - ref_pts_[i - 1];
+ Eigen::Vector3d v = (ref_pts_[i] - ref_pts_[i - 2]).normalized();
+
+ double xdotv = x.dot(v);
+ Eigen::Vector3d xcrossv = x.cross(v);
+
+ double f = pow((xdotv), 2) / a2 + pow(xcrossv.norm(), 2) / b2;
+ cost += f;
+
+ Eigen::Matrix3d m;
+ m << 0, -v(2), v(1), v(2), 0, -v(0), -v(1), v(0), 0;
+ Eigen::Vector3d df_dx = 2 * xdotv / a2 * v + 2 / b2 * m * xcrossv;
+
+ gradient.col(i - 1) += df_dx / 6;
+ gradient.col(i) += 4 * df_dx / 6;
+ gradient.col(i + 1) += df_dx / 6;
+ }
+ }
+
+ void BsplineOptimizer::calcSmoothnessCost(const Eigen::MatrixXd &q, double &cost,
+ Eigen::MatrixXd &gradient, bool falg_use_jerk /* = true*/)
+ {
+
+ cost = 0.0;
+
+ if (falg_use_jerk)
+ {
+ Eigen::Vector3d jerk, temp_j;
+
+ for (int i = 0; i < q.cols() - 3; i++)
+ {
+ /* evaluate jerk */
+ jerk = q.col(i + 3) - 3 * q.col(i + 2) + 3 * q.col(i + 1) - q.col(i);
+ cost += jerk.squaredNorm();
+ temp_j = 2.0 * jerk;
+ /* jerk gradient */
+ gradient.col(i + 0) += -temp_j;
+ gradient.col(i + 1) += 3.0 * temp_j;
+ gradient.col(i + 2) += -3.0 * temp_j;
+ gradient.col(i + 3) += temp_j;
+ }
+ }
+ else
+ {
+ Eigen::Vector3d acc, temp_acc;
+
+ for (int i = 0; i < q.cols() - 2; i++)
+ {
+ /* evaluate acc */
+ acc = q.col(i + 2) - 2 * q.col(i + 1) + q.col(i);
+ cost += acc.squaredNorm();
+ temp_acc = 2.0 * acc;
+ /* acc gradient */
+ gradient.col(i + 0) += temp_acc;
+ gradient.col(i + 1) += -2.0 * temp_acc;
+ gradient.col(i + 2) += temp_acc;
+ }
+ }
+ }
+
+ void BsplineOptimizer::calcTerminalCost(const Eigen::MatrixXd &q, double &cost, Eigen::MatrixXd &gradient)
+ {
+ cost = 0.0;
+
+ // zero cost and gradient in hard constraints
+ Eigen::Vector3d q_3, q_2, q_1, dq;
+ q_3 = q.col(q.cols()-3);
+ q_2 = q.col(q.cols()-2);
+ q_1 = q.col(q.cols()-1);
+
+ dq = 1 / 6.0 * (q_3 + 4 * q_2 + q_1) - local_target_pt_;
+ cost += dq.squaredNorm();
+
+ gradient.col(q.cols()-3) += 2 * dq * (1 / 6.0);
+ gradient.col(q.cols()-2) += 2 * dq * (4 / 6.0);
+ gradient.col(q.cols()-1) += 2 * dq * (1 / 6.0);
+
+ }
+
+ void BsplineOptimizer::calcFeasibilityCost(const Eigen::MatrixXd &q, double &cost,
+ Eigen::MatrixXd &gradient)
+ {
+
+ //#define SECOND_DERIVATIVE_CONTINOUS
+
+#ifdef SECOND_DERIVATIVE_CONTINOUS
+
+ cost = 0.0;
+ double demarcation = 1.0; // 1m/s, 1m/s/s
+ double ar = 3 * demarcation, br = -3 * pow(demarcation, 2), cr = pow(demarcation, 3);
+ double al = ar, bl = -br, cl = cr;
+
+ /* abbreviation */
+ double ts, ts_inv2, ts_inv3;
+ ts = bspline_interval_;
+ ts_inv2 = 1 / ts / ts;
+ ts_inv3 = 1 / ts / ts / ts;
+
+ /* velocity feasibility */
+ for (int i = 0; i < q.cols() - 1; i++)
+ {
+ Eigen::Vector3d vi = (q.col(i + 1) - q.col(i)) / ts;
+
+ for (int j = 0; j < 3; j++)
+ {
+ if (vi(j) > max_vel_ + demarcation)
+ {
+ double diff = vi(j) - max_vel_;
+ cost += (ar * diff * diff + br * diff + cr) * ts_inv3; // multiply ts_inv3 to make vel and acc has similar magnitude
+
+ double grad = (2.0 * ar * diff + br) / ts * ts_inv3;
+ gradient(j, i + 0) += -grad;
+ gradient(j, i + 1) += grad;
+ }
+ else if (vi(j) > max_vel_)
+ {
+ double diff = vi(j) - max_vel_;
+ cost += pow(diff, 3) * ts_inv3;
+ ;
+
+ double grad = 3 * diff * diff / ts * ts_inv3;
+ ;
+ gradient(j, i + 0) += -grad;
+ gradient(j, i + 1) += grad;
+ }
+ else if (vi(j) < -(max_vel_ + demarcation))
+ {
+ double diff = vi(j) + max_vel_;
+ cost += (al * diff * diff + bl * diff + cl) * ts_inv3;
+
+ double grad = (2.0 * al * diff + bl) / ts * ts_inv3;
+ gradient(j, i + 0) += -grad;
+ gradient(j, i + 1) += grad;
+ }
+ else if (vi(j) < -max_vel_)
+ {
+ double diff = vi(j) + max_vel_;
+ cost += -pow(diff, 3) * ts_inv3;
+
+ double grad = -3 * diff * diff / ts * ts_inv3;
+ gradient(j, i + 0) += -grad;
+ gradient(j, i + 1) += grad;
+ }
+ else
+ {
+ /* nothing happened */
+ }
+ }
+ }
+
+ /* acceleration feasibility */
+ for (int i = 0; i < q.cols() - 2; i++)
+ {
+ Eigen::Vector3d ai = (q.col(i + 2) - 2 * q.col(i + 1) + q.col(i)) * ts_inv2;
+
+ for (int j = 0; j < 3; j++)
+ {
+ if (ai(j) > max_acc_ + demarcation)
+ {
+ double diff = ai(j) - max_acc_;
+ cost += ar * diff * diff + br * diff + cr;
+
+ double grad = (2.0 * ar * diff + br) * ts_inv2;
+ gradient(j, i + 0) += grad;
+ gradient(j, i + 1) += -2 * grad;
+ gradient(j, i + 2) += grad;
+ }
+ else if (ai(j) > max_acc_)
+ {
+ double diff = ai(j) - max_acc_;
+ cost += pow(diff, 3);
+
+ double grad = 3 * diff * diff * ts_inv2;
+ gradient(j, i + 0) += grad;
+ gradient(j, i + 1) += -2 * grad;
+ gradient(j, i + 2) += grad;
+ }
+ else if (ai(j) < -(max_acc_ + demarcation))
+ {
+ double diff = ai(j) + max_acc_;
+ cost += al * diff * diff + bl * diff + cl;
+
+ double grad = (2.0 * al * diff + bl) * ts_inv2;
+ gradient(j, i + 0) += grad;
+ gradient(j, i + 1) += -2 * grad;
+ gradient(j, i + 2) += grad;
+ }
+ else if (ai(j) < -max_acc_)
+ {
+ double diff = ai(j) + max_acc_;
+ cost += -pow(diff, 3);
+
+ double grad = -3 * diff * diff * ts_inv2;
+ gradient(j, i + 0) += grad;
+ gradient(j, i + 1) += -2 * grad;
+ gradient(j, i + 2) += grad;
+ }
+ else
+ {
+ /* nothing happened */
+ }
+ }
+ }
+
+#else
+
+ cost = 0.0;
+ /* abbreviation */
+ double ts, /*vm2, am2, */ ts_inv2;
+ // vm2 = max_vel_ * max_vel_;
+ // am2 = max_acc_ * max_acc_;
+
+ ts = bspline_interval_;
+ ts_inv2 = 1 / ts / ts;
+
+ /* velocity feasibility */
+ for (int i = 0; i < q.cols() - 1; i++)
+ {
+ Eigen::Vector3d vi = (q.col(i + 1) - q.col(i)) / ts;
+
+ //cout << "temp_v * vi=" ;
+ for (int j = 0; j < 3; j++)
+ {
+ if (vi(j) > max_vel_)
+ {
+ // cout << "zx-todo VEL" << endl;
+ // cout << vi(j) << endl;
+ cost += pow(vi(j) - max_vel_, 2) * ts_inv2; // multiply ts_inv3 to make vel and acc has similar magnitude
+
+ gradient(j, i + 0) += -2 * (vi(j) - max_vel_) / ts * ts_inv2;
+ gradient(j, i + 1) += 2 * (vi(j) - max_vel_) / ts * ts_inv2;
+ }
+ else if (vi(j) < -max_vel_)
+ {
+ cost += pow(vi(j) + max_vel_, 2) * ts_inv2;
+
+ gradient(j, i + 0) += -2 * (vi(j) + max_vel_) / ts * ts_inv2;
+ gradient(j, i + 1) += 2 * (vi(j) + max_vel_) / ts * ts_inv2;
+ }
+ else
+ {
+ /* code */
+ }
+ }
+ }
+
+ /* acceleration feasibility */
+ for (int i = 0; i < q.cols() - 2; i++)
+ {
+ Eigen::Vector3d ai = (q.col(i + 2) - 2 * q.col(i + 1) + q.col(i)) * ts_inv2;
+
+ //cout << "temp_a * ai=" ;
+ for (int j = 0; j < 3; j++)
+ {
+ if (ai(j) > max_acc_)
+ {
+ // cout << "zx-todo ACC" << endl;
+ // cout << ai(j) << endl;
+ cost += pow(ai(j) - max_acc_, 2);
+
+ gradient(j, i + 0) += 2 * (ai(j) - max_acc_) * ts_inv2;
+ gradient(j, i + 1) += -4 * (ai(j) - max_acc_) * ts_inv2;
+ gradient(j, i + 2) += 2 * (ai(j) - max_acc_) * ts_inv2;
+ }
+ else if (ai(j) < -max_acc_)
+ {
+ cost += pow(ai(j) + max_acc_, 2);
+
+ gradient(j, i + 0) += 2 * (ai(j) + max_acc_) * ts_inv2;
+ gradient(j, i + 1) += -4 * (ai(j) + max_acc_) * ts_inv2;
+ gradient(j, i + 2) += 2 * (ai(j) + max_acc_) * ts_inv2;
+ }
+ else
+ {
+ /* code */
+ }
+ }
+ //cout << endl;
+ }
+
+#endif
+ }
+
+ bool BsplineOptimizer::check_collision_and_rebound(void)
+ {
+
+ int end_idx = cps_.size - order_;
+
+ /*** Check and segment the initial trajectory according to obstacles ***/
+ int in_id, out_id;
+ vector> segment_ids;
+ bool flag_new_obs_valid = false;
+ int i_end = end_idx - (end_idx - order_) / 3;
+ for (int i = order_ - 1; i <= i_end; ++i)
+ {
+
+ bool occ = grid_map_->getInflateOccupancy(cps_.points.col(i));
+
+ /*** check if the new collision will be valid ***/
+ if (occ)
+ {
+ for (size_t k = 0; k < cps_.direction[i].size(); ++k)
+ {
+ cout.precision(2);
+ if ((cps_.points.col(i) - cps_.base_point[i][k]).dot(cps_.direction[i][k]) < 1 * grid_map_->getResolution()) // current point is outside all the collision_points.
+ {
+ occ = false; // Not really takes effect, just for better hunman understanding.
+ break;
+ }
+ }
+ }
+
+ if (occ)
+ {
+ flag_new_obs_valid = true;
+
+ int j;
+ for (j = i - 1; j >= 0; --j)
+ {
+ occ = grid_map_->getInflateOccupancy(cps_.points.col(j));
+ if (!occ)
+ {
+ in_id = j;
+ break;
+ }
+ }
+ if (j < 0) // fail to get the obs free point
+ {
+ ROS_ERROR("ERROR! the drone is in obstacle. This should not happen.");
+ in_id = 0;
+ }
+
+ for (j = i + 1; j < cps_.size; ++j)
+ {
+ occ = grid_map_->getInflateOccupancy(cps_.points.col(j));
+
+ if (!occ)
+ {
+ out_id = j;
+ break;
+ }
+ }
+ if (j >= cps_.size) // fail to get the obs free point
+ {
+ ROS_WARN("WARN! terminal point of the current trajectory is in obstacle, skip this planning.");
+
+ force_stop_type_ = STOP_FOR_ERROR;
+ return false;
+ }
+
+ i = j + 1;
+
+ segment_ids.push_back(std::pair(in_id, out_id));
+ }
+ }
+
+ if (flag_new_obs_valid)
+ {
+ vector> a_star_pathes;
+ for (size_t i = 0; i < segment_ids.size(); ++i)
+ {
+ /*** a star search ***/
+ Eigen::Vector3d in(cps_.points.col(segment_ids[i].first)), out(cps_.points.col(segment_ids[i].second));
+ if (a_star_->AstarSearch(/*(in-out).norm()/10+0.05*/ 0.1, in, out))
+ {
+ a_star_pathes.push_back(a_star_->getPath());
+ }
+ else
+ {
+ ROS_ERROR("a star error");
+ segment_ids.erase(segment_ids.begin() + i);
+ i--;
+ }
+ }
+
+ for (size_t i = 1; i < segment_ids.size(); i++) // Avoid overlap
+ {
+ if (segment_ids[i - 1].second >= segment_ids[i].first)
+ {
+ double middle = (double)(segment_ids[i - 1].second + segment_ids[i].first) / 2.0;
+ segment_ids[i - 1].second = static_cast(middle - 0.1);
+ segment_ids[i].first = static_cast(middle + 1.1);
+ }
+ }
+
+ /*** Assign parameters to each segment ***/
+ for (size_t i = 0; i < segment_ids.size(); ++i)
+ {
+ // step 1
+ for (int j = segment_ids[i].first; j <= segment_ids[i].second; ++j)
+ cps_.flag_temp[j] = false;
+
+ // step 2
+ int got_intersection_id = -1;
+ for (int j = segment_ids[i].first + 1; j < segment_ids[i].second; ++j)
+ {
+ Eigen::Vector3d ctrl_pts_law(cps_.points.col(j + 1) - cps_.points.col(j - 1)), intersection_point;
+ int Astar_id = a_star_pathes[i].size() / 2, last_Astar_id; // Let "Astar_id = id_of_the_most_far_away_Astar_point" will be better, but it needs more computation
+ double val = (a_star_pathes[i][Astar_id] - cps_.points.col(j)).dot(ctrl_pts_law), last_val = val;
+ while (Astar_id >= 0 && Astar_id < (int)a_star_pathes[i].size())
+ {
+ last_Astar_id = Astar_id;
+
+ if (val >= 0)
+ --Astar_id;
+ else
+ ++Astar_id;
+
+ val = (a_star_pathes[i][Astar_id] - cps_.points.col(j)).dot(ctrl_pts_law);
+
+ // cout << val << endl;
+
+ if (val * last_val <= 0 && (abs(val) > 0 || abs(last_val) > 0)) // val = last_val = 0.0 is not allowed
+ {
+ intersection_point =
+ a_star_pathes[i][Astar_id] +
+ ((a_star_pathes[i][Astar_id] - a_star_pathes[i][last_Astar_id]) *
+ (ctrl_pts_law.dot(cps_.points.col(j) - a_star_pathes[i][Astar_id]) / ctrl_pts_law.dot(a_star_pathes[i][Astar_id] - a_star_pathes[i][last_Astar_id])) // = t
+ );
+
+ got_intersection_id = j;
+ break;
+ }
+ }
+
+ if (got_intersection_id >= 0)
+ {
+ double length = (intersection_point - cps_.points.col(j)).norm();
+ if (length > 1e-5)
+ {
+ cps_.flag_temp[j] = true;
+ for (double a = length; a >= 0.0; a -= grid_map_->getResolution())
+ {
+ bool occ = grid_map_->getInflateOccupancy((a / length) * intersection_point + (1 - a / length) * cps_.points.col(j));
+
+ if (occ || a < grid_map_->getResolution())
+ {
+ if (occ)
+ a += grid_map_->getResolution();
+ cps_.base_point[j].push_back((a / length) * intersection_point + (1 - a / length) * cps_.points.col(j));
+ cps_.direction[j].push_back((intersection_point - cps_.points.col(j)).normalized());
+ break;
+ }
+ }
+ }
+ else
+ {
+ got_intersection_id = -1;
+ }
+ }
+ }
+
+ //step 3
+ if (got_intersection_id >= 0)
+ {
+ for (int j = got_intersection_id + 1; j <= segment_ids[i].second; ++j)
+ if (!cps_.flag_temp[j])
+ {
+ cps_.base_point[j].push_back(cps_.base_point[j - 1].back());
+ cps_.direction[j].push_back(cps_.direction[j - 1].back());
+ }
+
+ for (int j = got_intersection_id - 1; j >= segment_ids[i].first; --j)
+ if (!cps_.flag_temp[j])
+ {
+ cps_.base_point[j].push_back(cps_.base_point[j + 1].back());
+ cps_.direction[j].push_back(cps_.direction[j + 1].back());
+ }
+ }
+ else
+ ROS_WARN("Failed to generate direction. It doesn't matter.");
+ }
+
+ force_stop_type_ = STOP_FOR_REBOUND;
+ return true;
+ }
+
+ return false;
+ }
+
+ bool BsplineOptimizer::BsplineOptimizeTrajRebound(Eigen::MatrixXd &optimal_points, double ts)
+ {
+ setBsplineInterval(ts);
+
+ double final_cost;
+ bool flag_success = rebound_optimize(final_cost);
+
+ optimal_points = cps_.points;
+
+ return flag_success;
+ }
+
+ bool BsplineOptimizer::BsplineOptimizeTrajRebound(Eigen::MatrixXd &optimal_points, double &final_cost, const ControlPoints &control_points, double ts)
+ {
+ setBsplineInterval(ts);
+
+ cps_ = control_points;
+
+ bool flag_success = rebound_optimize(final_cost);
+
+ optimal_points = cps_.points;
+
+ return flag_success;
+ }
+
+ bool BsplineOptimizer::BsplineOptimizeTrajRefine(const Eigen::MatrixXd &init_points, const double ts, Eigen::MatrixXd &optimal_points)
+ {
+
+ setControlPoints(init_points);
+ setBsplineInterval(ts);
+
+ bool flag_success = refine_optimize();
+
+ optimal_points = cps_.points;
+
+ return flag_success;
+ }
+
+ bool BsplineOptimizer::rebound_optimize(double &final_cost)
+ {
+ iter_num_ = 0;
+ int start_id = order_;
+ // int end_id = this->cps_.size - order_; //Fixed end
+ int end_id = this->cps_.size; // Free end
+ variable_num_ = 3 * (end_id - start_id);
+
+ ros::Time t0 = ros::Time::now(), t1, t2;
+ int restart_nums = 0, rebound_times = 0;
+ ;
+ bool flag_force_return, flag_occ, success;
+ new_lambda2_ = lambda2_;
+ constexpr int MAX_RESART_NUMS_SET = 3;
+ do
+ {
+ /* ---------- prepare ---------- */
+ min_cost_ = std::numeric_limits::max();
+ min_ellip_dist_ = INIT_min_ellip_dist_;
+ iter_num_ = 0;
+ flag_force_return = false;
+ flag_occ = false;
+ success = false;
+
+ double q[variable_num_];
+ memcpy(q, cps_.points.data() + 3 * start_id, variable_num_ * sizeof(q[0]));
+
+ lbfgs::lbfgs_parameter_t lbfgs_params;
+ lbfgs::lbfgs_load_default_parameters(&lbfgs_params);
+ lbfgs_params.mem_size = 16;
+ lbfgs_params.max_iterations = 200;
+ lbfgs_params.g_epsilon = 0.01;
+
+ /* ---------- optimize ---------- */
+ t1 = ros::Time::now();
+ int result = lbfgs::lbfgs_optimize(variable_num_, q, &final_cost, BsplineOptimizer::costFunctionRebound, NULL, BsplineOptimizer::earlyExit, this, &lbfgs_params);
+ t2 = ros::Time::now();
+ double time_ms = (t2 - t1).toSec() * 1000;
+ double total_time_ms = (t2 - t0).toSec() * 1000;
+
+ /* ---------- success temporary, check collision again ---------- */
+ if (result == lbfgs::LBFGS_CONVERGENCE ||
+ result == lbfgs::LBFGSERR_MAXIMUMITERATION ||
+ result == lbfgs::LBFGS_ALREADY_MINIMIZED ||
+ result == lbfgs::LBFGS_STOP)
+ {
+ //ROS_WARN("Solver error in planning!, return = %s", lbfgs::lbfgs_strerror(result));
+ flag_force_return = false;
+
+ /*** collision check, phase 1 ***/
+ if ((min_ellip_dist_ != INIT_min_ellip_dist_) && (min_ellip_dist_ > swarm_clearance_))
+ {
+ success = false;
+ restart_nums++;
+ initControlPoints(cps_.points, false);
+ new_lambda2_ *= 2;
+
+ printf("\033[32miter(+1)=%d,time(ms)=%5.3f, swarm too close, keep optimizing\n\033[0m", iter_num_, time_ms);
+
+ continue;
+ }
+
+ /*** collision check, phase 2 ***/
+ UniformBspline traj = UniformBspline(cps_.points, 3, bspline_interval_);
+ double tm, tmp;
+ traj.getTimeSpan(tm, tmp);
+ double t_step = (tmp - tm) / ((traj.evaluateDeBoorT(tmp) - traj.evaluateDeBoorT(tm)).norm() / grid_map_->getResolution());
+ for (double t = tm; t < tmp * 2 / 3; t += t_step) // Only check the closest 2/3 partition of the whole trajectory.
+ {
+ flag_occ = grid_map_->getInflateOccupancy(traj.evaluateDeBoorT(t));
+ if (flag_occ)
+ {
+ //cout << "hit_obs, t=" << t << " P=" << traj.evaluateDeBoorT(t).transpose() << endl;
+
+ if (t <= bspline_interval_) // First 3 control points in obstacles!
+ {
+ // cout << cps_.points.col(1).transpose() << "\n"
+ // << cps_.points.col(2).transpose() << "\n"
+ // << cps_.points.col(3).transpose() << "\n"
+ // << cps_.points.col(4).transpose() << endl;
+ ROS_WARN("First 3 control points in obstacles! return false, t=%f", t);
+ return false;
+ }
+
+ break;
+ }
+ }
+
+ // cout << "XXXXXX" << ((cps_.points.col(cps_.points.cols()-1) + 4*cps_.points.col(cps_.points.cols()-2) + cps_.points.col(cps_.points.cols()-3))/6 - local_target_pt_).norm() << endl;
+
+ /*** collision check, phase 3 ***/
+//#define USE_SECOND_CLEARENCE_CHECK
+#ifdef USE_SECOND_CLEARENCE_CHECK
+ bool flag_cls_xyp, flag_cls_xyn, flag_cls_zp, flag_cls_zn;
+ Eigen::Vector3d start_end_vec = traj.evaluateDeBoorT(tmp) - traj.evaluateDeBoorT(tm);
+ Eigen::Vector3d offset_xy(-start_end_vec(0), start_end_vec(1), 0);
+ offset_xy.normalize();
+ Eigen::Vector3d offset_z = start_end_vec.cross(offset_xy);
+ offset_z.normalize();
+ offset_xy *= cps_.clearance / 2;
+ offset_z *= cps_.clearance / 2;
+
+ Eigen::MatrixXd check_pts(cps_.points.rows(), cps_.points.cols());
+ for (Eigen::Index i = 0; i < cps_.points.cols(); i++)
+ {
+ check_pts.col(i) = cps_.points.col(i);
+ check_pts(0, i) += offset_xy(0);
+ check_pts(1, i) += offset_xy(1);
+ check_pts(2, i) += offset_xy(2);
+ }
+ flag_cls_xyp = initControlPoints(check_pts, false).size() > 0;
+ for (Eigen::Index i = 0; i < cps_.points.cols(); i++)
+ {
+ check_pts(0, i) -= 2 * offset_xy(0);
+ check_pts(1, i) -= 2 * offset_xy(1);
+ check_pts(2, i) -= 2 * offset_xy(2);
+ }
+ flag_cls_xyn = initControlPoints(check_pts, false).size() > 0;
+ for (Eigen::Index i = 0; i < cps_.points.cols(); i++)
+ {
+ check_pts(0, i) += offset_xy(0) + offset_z(0);
+ check_pts(1, i) += offset_xy(1) + offset_z(1);
+ check_pts(2, i) += offset_xy(2) + offset_z(2);
+ }
+ flag_cls_zp = initControlPoints(check_pts, false).size() > 0;
+ for (Eigen::Index i = 0; i < cps_.points.cols(); i++)
+ {
+ check_pts(0, i) -= 2 * offset_z(0);
+ check_pts(1, i) -= 2 * offset_z(1);
+ check_pts(2, i) -= 2 * offset_z(2);
+ }
+ flag_cls_zn = initControlPoints(check_pts, false).size() > 0;
+ if ((flag_cls_xyp ^ flag_cls_xyn) || (flag_cls_zp ^ flag_cls_zn))
+ flag_occ = true;
+#endif
+
+ if (!flag_occ)
+ {
+ printf("\033[32miter(+1)=%d,time(ms)=%5.3f,total_t(ms)=%5.3f,cost=%5.3f\n\033[0m", iter_num_, time_ms, total_time_ms, final_cost);
+ success = true;
+ }
+ else // restart
+ {
+ restart_nums++;
+ initControlPoints(cps_.points, false);
+ new_lambda2_ *= 2;
+
+ printf("\033[32miter(+1)=%d,time(ms)=%5.3f, collided, keep optimizing\n\033[0m", iter_num_, time_ms);
+ }
+ }
+ else if (result == lbfgs::LBFGSERR_CANCELED)
+ {
+ flag_force_return = true;
+ rebound_times++;
+ cout << "iter=" << iter_num_ << ",time(ms)=" << time_ms << ",rebound." << endl;
+ }
+ else
+ {
+ ROS_WARN("Solver error. Return = %d, %s. Skip this planning.", result, lbfgs::lbfgs_strerror(result));
+ // while (ros::ok());
+ }
+
+ } while (
+ ((flag_occ || ((min_ellip_dist_ != INIT_min_ellip_dist_) && (min_ellip_dist_ > swarm_clearance_))) && restart_nums < MAX_RESART_NUMS_SET) ||
+ (flag_force_return && force_stop_type_ == STOP_FOR_REBOUND && rebound_times <= 20));
+
+ return success;
+ }
+
+ bool BsplineOptimizer::refine_optimize()
+ {
+ iter_num_ = 0;
+ int start_id = order_;
+ int end_id = this->cps_.points.cols() - order_;
+ variable_num_ = 3 * (end_id - start_id);
+
+ double q[variable_num_];
+ double final_cost;
+
+ memcpy(q, cps_.points.data() + 3 * start_id, variable_num_ * sizeof(q[0]));
+
+ double origin_lambda4 = lambda4_;
+ bool flag_safe = true;
+ int iter_count = 0;
+ do
+ {
+ lbfgs::lbfgs_parameter_t lbfgs_params;
+ lbfgs::lbfgs_load_default_parameters(&lbfgs_params);
+ lbfgs_params.mem_size = 16;
+ lbfgs_params.max_iterations = 200;
+ lbfgs_params.g_epsilon = 0.001;
+
+ int result = lbfgs::lbfgs_optimize(variable_num_, q, &final_cost, BsplineOptimizer::costFunctionRefine, NULL, NULL, this, &lbfgs_params);
+ if (result == lbfgs::LBFGS_CONVERGENCE ||
+ result == lbfgs::LBFGSERR_MAXIMUMITERATION ||
+ result == lbfgs::LBFGS_ALREADY_MINIMIZED ||
+ result == lbfgs::LBFGS_STOP)
+ {
+ //pass
+ }
+ else
+ {
+ ROS_ERROR("Solver error in refining!, return = %d, %s", result, lbfgs::lbfgs_strerror(result));
+ }
+
+ UniformBspline traj = UniformBspline(cps_.points, 3, bspline_interval_);
+ double tm, tmp;
+ traj.getTimeSpan(tm, tmp);
+ double t_step = (tmp - tm) / ((traj.evaluateDeBoorT(tmp) - traj.evaluateDeBoorT(tm)).norm() / grid_map_->getResolution()); // Step size is defined as the maximum size that can passes throgth every gird.
+ for (double t = tm; t < tmp * 2 / 3; t += t_step)
+ {
+ if (grid_map_->getInflateOccupancy(traj.evaluateDeBoorT(t)))
+ {
+ // cout << "Refined traj hit_obs, t=" << t << " P=" << traj.evaluateDeBoorT(t).transpose() << endl;
+
+ Eigen::MatrixXd ref_pts(ref_pts_.size(), 3);
+ for (size_t i = 0; i < ref_pts_.size(); i++)
+ {
+ ref_pts.row(i) = ref_pts_[i].transpose();
+ }
+
+ flag_safe = false;
+ break;
+ }
+ }
+
+ if (!flag_safe)
+ lambda4_ *= 2;
+
+ iter_count++;
+ } while (!flag_safe && iter_count <= 0);
+
+ lambda4_ = origin_lambda4;
+
+ //cout << "iter_num_=" << iter_num_ << endl;
+
+ return flag_safe;
+ }
+
+ void BsplineOptimizer::combineCostRebound(const double *x, double *grad, double &f_combine, const int n)
+ {
+ // cout << "drone_id_=" << drone_id_ << endl;
+ // cout << "cps_.points.size()=" << cps_.points.size() << endl;
+ // cout << "n=" << n << endl;
+ // cout << "sizeof(x[0])=" << sizeof(x[0]) << endl;
+
+ memcpy(cps_.points.data() + 3 * order_, x, n * sizeof(x[0]));
+
+ /* ---------- evaluate cost and gradient ---------- */
+ double f_smoothness, f_distance, f_feasibility /*, f_mov_objs*/, f_swarm, f_terminal;
+
+ Eigen::MatrixXd g_smoothness = Eigen::MatrixXd::Zero(3, cps_.size);
+ Eigen::MatrixXd g_distance = Eigen::MatrixXd::Zero(3, cps_.size);
+ Eigen::MatrixXd g_feasibility = Eigen::MatrixXd::Zero(3, cps_.size);
+ // Eigen::MatrixXd g_mov_objs = Eigen::MatrixXd::Zero(3, cps_.size);
+ Eigen::MatrixXd g_swarm = Eigen::MatrixXd::Zero(3, cps_.size);
+ Eigen::MatrixXd g_terminal = Eigen::MatrixXd::Zero(3, cps_.size);
+
+ calcSmoothnessCost(cps_.points, f_smoothness, g_smoothness);
+ calcDistanceCostRebound(cps_.points, f_distance, g_distance, iter_num_, f_smoothness);
+ calcFeasibilityCost(cps_.points, f_feasibility, g_feasibility);
+ // calcMovingObjCost(cps_.points, f_mov_objs, g_mov_objs);
+ calcSwarmCost(cps_.points, f_swarm, g_swarm);
+ calcTerminalCost( cps_.points, f_terminal, g_terminal );
+
+ f_combine = lambda1_ * f_smoothness + new_lambda2_ * f_distance + lambda3_ * f_feasibility + new_lambda2_ * f_swarm + lambda2_ * f_terminal;
+ //f_combine = lambda1_ * f_smoothness + new_lambda2_ * f_distance + lambda3_ * f_feasibility + new_lambda2_ * f_mov_objs;
+ //printf("origin %f %f %f %f\n", f_smoothness, f_distance, f_feasibility, f_combine);
+
+ Eigen::MatrixXd grad_3D = lambda1_ * g_smoothness + new_lambda2_ * g_distance + lambda3_ * g_feasibility + new_lambda2_ * g_swarm + lambda2_ * g_terminal;
+ //Eigen::MatrixXd grad_3D = lambda1_ * g_smoothness + new_lambda2_ * g_distance + lambda3_ * g_feasibility + new_lambda2_ * g_mov_objs;
+ memcpy(grad, grad_3D.data() + 3 * order_, n * sizeof(grad[0]));
+ }
+
+ void BsplineOptimizer::combineCostRefine(const double *x, double *grad, double &f_combine, const int n)
+ {
+
+ memcpy(cps_.points.data() + 3 * order_, x, n * sizeof(x[0]));
+
+ /* ---------- evaluate cost and gradient ---------- */
+ double f_smoothness, f_fitness, f_feasibility;
+
+ Eigen::MatrixXd g_smoothness = Eigen::MatrixXd::Zero(3, cps_.points.cols());
+ Eigen::MatrixXd g_fitness = Eigen::MatrixXd::Zero(3, cps_.points.cols());
+ Eigen::MatrixXd g_feasibility = Eigen::MatrixXd::Zero(3, cps_.points.cols());
+
+ //time_satrt = ros::Time::now();
+
+ calcSmoothnessCost(cps_.points, f_smoothness, g_smoothness);
+ calcFitnessCost(cps_.points, f_fitness, g_fitness);
+ calcFeasibilityCost(cps_.points, f_feasibility, g_feasibility);
+
+ /* ---------- convert to solver format...---------- */
+ f_combine = lambda1_ * f_smoothness + lambda4_ * f_fitness + lambda3_ * f_feasibility;
+ // printf("origin %f %f %f %f\n", f_smoothness, f_fitness, f_feasibility, f_combine);
+
+ Eigen::MatrixXd grad_3D = lambda1_ * g_smoothness + lambda4_ * g_fitness + lambda3_ * g_feasibility;
+ memcpy(grad, grad_3D.data() + 3 * order_, n * sizeof(grad[0]));
+ }
+
+} // namespace ego_planner
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/bspline_opt/src/gradient_descent_optimizer.cpp b/motion_planning/3d/ego_planner/planner/bspline_opt/src/gradient_descent_optimizer.cpp
new file mode 100755
index 0000000..c3371d6
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/bspline_opt/src/gradient_descent_optimizer.cpp
@@ -0,0 +1,94 @@
+#include
+
+#define RESET "\033[0m"
+#define RED "\033[31m"
+
+GradientDescentOptimizer::RESULT
+GradientDescentOptimizer::optimize(Eigen::VectorXd &x_init_optimal, double &opt_f)
+{
+ if (min_grad_ < 1e-10)
+ {
+ cout << RED << "min_grad_ is invalid:" << min_grad_ << RESET << endl;
+ return FAILED;
+ }
+ if (iter_limit_ <= 2)
+ {
+ cout << RED << "iter_limit_ is invalid:" << iter_limit_ << RESET << endl;
+ return FAILED;
+ }
+
+ void *f_data = f_data_;
+ int iter = 2;
+ int invoke_count = 2;
+ bool force_return;
+ Eigen::VectorXd x_k(x_init_optimal), x_kp1(x_init_optimal.rows());
+ double cost_k, cost_kp1, cost_min;
+ Eigen::VectorXd grad_k(x_init_optimal.rows()), grad_kp1(x_init_optimal.rows());
+
+ cost_k = objfun_(x_k, grad_k, force_return, f_data);
+ if (force_return)
+ return RETURN_BY_ORDER;
+ cost_min = cost_k;
+ double max_grad = max(abs(grad_k.maxCoeff()), abs(grad_k.minCoeff()));
+ constexpr double MAX_MOVEMENT_AT_FIRST_ITERATION = 0.1; // meter
+ double alpha0 = max_grad < MAX_MOVEMENT_AT_FIRST_ITERATION ? 1.0 : (MAX_MOVEMENT_AT_FIRST_ITERATION / max_grad);
+ x_kp1 = x_k - alpha0 * grad_k;
+ cost_kp1 = objfun_(x_kp1, grad_kp1, force_return, f_data);
+ if (force_return)
+ return RETURN_BY_ORDER;
+ if (cost_min > cost_kp1)
+ cost_min = cost_kp1;
+
+ /*** start iteration ***/
+ while (++iter <= iter_limit_ && invoke_count <= invoke_limit_)
+ {
+ Eigen::VectorXd s = x_kp1 - x_k;
+ Eigen::VectorXd y = grad_kp1 - grad_k;
+ double alpha = s.dot(y) / y.dot(y);
+ if (isnan(alpha) || isinf(alpha))
+ {
+ cout << RED << "step size invalid! alpha=" << alpha << RESET << endl;
+ return FAILED;
+ }
+
+ if (iter % 2) // to aviod copying operations
+ {
+ do
+ {
+ x_k = x_kp1 - alpha * grad_kp1;
+ cost_k = objfun_(x_k, grad_k, force_return, f_data);
+ invoke_count++;
+ if (force_return)
+ return RETURN_BY_ORDER;
+ alpha *= 0.5;
+ } while (cost_k > cost_kp1 - 1e-4 * alpha * grad_kp1.transpose() * grad_kp1); // Armijo condition
+
+ if (grad_k.norm() < min_grad_)
+ {
+ opt_f = cost_k;
+ return FIND_MIN;
+ }
+ }
+ else
+ {
+ do
+ {
+ x_kp1 = x_k - alpha * grad_k;
+ cost_kp1 = objfun_(x_kp1, grad_kp1, force_return, f_data);
+ invoke_count++;
+ if (force_return)
+ return RETURN_BY_ORDER;
+ alpha *= 0.5;
+ } while (cost_kp1 > cost_k - 1e-4 * alpha * grad_k.transpose() * grad_k); // Armijo condition
+
+ if (grad_kp1.norm() < min_grad_)
+ {
+ opt_f = cost_kp1;
+ return FIND_MIN;
+ }
+ }
+ }
+
+ opt_f = iter_limit_ % 2 ? cost_k : cost_kp1;
+ return REACH_MAX_ITERATION;
+}
diff --git a/motion_planning/3d/ego_planner/planner/bspline_opt/src/uniform_bspline.cpp b/motion_planning/3d/ego_planner/planner/bspline_opt/src/uniform_bspline.cpp
new file mode 100644
index 0000000..ff798f0
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/bspline_opt/src/uniform_bspline.cpp
@@ -0,0 +1,377 @@
+#include "bspline_opt/uniform_bspline.h"
+#include
+
+namespace ego_planner
+{
+
+ UniformBspline::UniformBspline(const Eigen::MatrixXd &points, const int &order,
+ const double &interval)
+ {
+ setUniformBspline(points, order, interval);
+ }
+
+ UniformBspline::~UniformBspline() {}
+
+ void UniformBspline::setUniformBspline(const Eigen::MatrixXd &points, const int &order,
+ const double &interval)
+ {
+ control_points_ = points;
+ p_ = order;
+ interval_ = interval;
+
+ n_ = points.cols() - 1;
+ m_ = n_ + p_ + 1;
+
+ u_ = Eigen::VectorXd::Zero(m_ + 1);
+ for (int i = 0; i <= m_; ++i)
+ {
+
+ if (i <= p_)
+ {
+ u_(i) = double(-p_ + i) * interval_;
+ }
+ else if (i > p_ && i <= m_ - p_)
+ {
+ u_(i) = u_(i - 1) + interval_;
+ }
+ else if (i > m_ - p_)
+ {
+ u_(i) = u_(i - 1) + interval_;
+ }
+ }
+ }
+
+ void UniformBspline::setKnot(const Eigen::VectorXd &knot) { this->u_ = knot; }
+
+ Eigen::VectorXd UniformBspline::getKnot() { return this->u_; }
+
+ bool UniformBspline::getTimeSpan(double &um, double &um_p)
+ {
+ if (p_ > u_.rows() || m_ - p_ > u_.rows())
+ return false;
+
+ um = u_(p_);
+ um_p = u_(m_ - p_);
+
+ return true;
+ }
+
+ Eigen::MatrixXd UniformBspline::getControlPoint() { return control_points_; }
+
+ Eigen::VectorXd UniformBspline::evaluateDeBoor(const double &u)
+ {
+
+ double ub = min(max(u_(p_), u), u_(m_ - p_));
+
+ // determine which [ui,ui+1] lay in
+ int k = p_;
+ while (true)
+ {
+ if (u_(k + 1) >= ub)
+ break;
+ ++k;
+ }
+
+ /* deBoor's alg */
+ vector d;
+ for (int i = 0; i <= p_; ++i)
+ {
+ d.push_back(control_points_.col(k - p_ + i));
+ // cout << d[i].transpose() << endl;
+ }
+
+ for (int r = 1; r <= p_; ++r)
+ {
+ for (int i = p_; i >= r; --i)
+ {
+ double alpha = (ub - u_[i + k - p_]) / (u_[i + 1 + k - r] - u_[i + k - p_]);
+ // cout << "alpha: " << alpha << endl;
+ d[i] = (1 - alpha) * d[i - 1] + alpha * d[i];
+ }
+ }
+
+ return d[p_];
+ }
+
+ // Eigen::VectorXd UniformBspline::evaluateDeBoorT(const double& t) {
+ // return evaluateDeBoor(t + u_(p_));
+ // }
+
+ Eigen::MatrixXd UniformBspline::getDerivativeControlPoints()
+ {
+ // The derivative of a b-spline is also a b-spline, its order become p_-1
+ // control point Qi = p_*(Pi+1-Pi)/(ui+p_+1-ui+1)
+ Eigen::MatrixXd ctp(control_points_.rows(), control_points_.cols() - 1);
+ for (int i = 0; i < ctp.cols(); ++i)
+ {
+ ctp.col(i) =
+ p_ * (control_points_.col(i + 1) - control_points_.col(i)) / (u_(i + p_ + 1) - u_(i + 1));
+ }
+ return ctp;
+ }
+
+ UniformBspline UniformBspline::getDerivative()
+ {
+ Eigen::MatrixXd ctp = getDerivativeControlPoints();
+ UniformBspline derivative(ctp, p_ - 1, interval_);
+
+ /* cut the first and last knot */
+ Eigen::VectorXd knot(u_.rows() - 2);
+ knot = u_.segment(1, u_.rows() - 2);
+ derivative.setKnot(knot);
+
+ return derivative;
+ }
+
+ double UniformBspline::getInterval() { return interval_; }
+
+ void UniformBspline::setPhysicalLimits(const double &vel, const double &acc, const double &tolerance)
+ {
+ limit_vel_ = vel;
+ limit_acc_ = acc;
+ limit_ratio_ = 1.1;
+ feasibility_tolerance_ = tolerance;
+ }
+
+ bool UniformBspline::checkFeasibility(double &ratio, bool show)
+ {
+ bool fea = true;
+
+ Eigen::MatrixXd P = control_points_;
+ int dimension = control_points_.rows();
+
+ /* check vel feasibility and insert points */
+ double max_vel = -1.0;
+ double enlarged_vel_lim = limit_vel_ * (1.0 + feasibility_tolerance_) + 1e-4;
+ for (int i = 0; i < P.cols() - 1; ++i)
+ {
+ Eigen::VectorXd vel = p_ * (P.col(i + 1) - P.col(i)) / (u_(i + p_ + 1) - u_(i + 1));
+
+ if (fabs(vel(0)) > enlarged_vel_lim || fabs(vel(1)) > enlarged_vel_lim ||
+ fabs(vel(2)) > enlarged_vel_lim)
+ {
+
+ if (show)
+ cout << "[Check]: Infeasible vel " << i << " :" << vel.transpose() << endl;
+ fea = false;
+
+ for (int j = 0; j < dimension; ++j)
+ {
+ max_vel = max(max_vel, fabs(vel(j)));
+ }
+ }
+ }
+
+ /* acc feasibility */
+ double max_acc = -1.0;
+ double enlarged_acc_lim = limit_acc_ * (1.0 + feasibility_tolerance_) + 1e-4;
+ for (int i = 0; i < P.cols() - 2; ++i)
+ {
+
+ Eigen::VectorXd acc = p_ * (p_ - 1) *
+ ((P.col(i + 2) - P.col(i + 1)) / (u_(i + p_ + 2) - u_(i + 2)) -
+ (P.col(i + 1) - P.col(i)) / (u_(i + p_ + 1) - u_(i + 1))) /
+ (u_(i + p_ + 1) - u_(i + 2));
+
+ if (fabs(acc(0)) > enlarged_acc_lim || fabs(acc(1)) > enlarged_acc_lim ||
+ fabs(acc(2)) > enlarged_acc_lim)
+ {
+
+ if (show)
+ cout << "[Check]: Infeasible acc " << i << " :" << acc.transpose() << endl;
+ fea = false;
+
+ for (int j = 0; j < dimension; ++j)
+ {
+ max_acc = max(max_acc, fabs(acc(j)));
+ }
+ }
+ }
+
+ ratio = max(max_vel / limit_vel_, sqrt(fabs(max_acc) / limit_acc_));
+
+ return fea;
+ }
+
+ void UniformBspline::lengthenTime(const double &ratio)
+ {
+ int num1 = 5;
+ int num2 = getKnot().rows() - 1 - 5;
+
+ double delta_t = (ratio - 1.0) * (u_(num2) - u_(num1));
+ double t_inc = delta_t / double(num2 - num1);
+ for (int i = num1 + 1; i <= num2; ++i)
+ u_(i) += double(i - num1) * t_inc;
+ for (int i = num2 + 1; i < u_.rows(); ++i)
+ u_(i) += delta_t;
+ }
+
+ // void UniformBspline::recomputeInit() {}
+
+ void UniformBspline::parameterizeToBspline(const double &ts, const vector &point_set,
+ const vector &start_end_derivative,
+ Eigen::MatrixXd &ctrl_pts)
+ {
+ if (ts <= 0)
+ {
+ cout << "[B-spline]:time step error." << endl;
+ return;
+ }
+
+ if (point_set.size() <= 3)
+ {
+ cout << "[B-spline]:point set have only " << point_set.size() << " points." << endl;
+ return;
+ }
+
+ if (start_end_derivative.size() != 4)
+ {
+ cout << "[B-spline]:derivatives error." << endl;
+ }
+
+ int K = point_set.size();
+
+ // write A
+ Eigen::Vector3d prow(3), vrow(3), arow(3);
+ prow << 1, 4, 1;
+ vrow << -1, 0, 1;
+ arow << 1, -2, 1;
+
+ Eigen::MatrixXd A = Eigen::MatrixXd::Zero(K + 4, K + 2);
+
+ for (int i = 0; i < K; ++i)
+ A.block(i, i, 1, 3) = (1 / 6.0) * prow.transpose();
+
+ A.block(K, 0, 1, 3) = (1 / 2.0 / ts) * vrow.transpose();
+ A.block(K + 1, K - 1, 1, 3) = (1 / 2.0 / ts) * vrow.transpose();
+
+ A.block(K + 2, 0, 1, 3) = (1 / ts / ts) * arow.transpose();
+ A.block(K + 3, K - 1, 1, 3) = (1 / ts / ts) * arow.transpose();
+
+ //cout << "A" << endl << A << endl << endl;
+
+ // write b
+ Eigen::VectorXd bx(K + 4), by(K + 4), bz(K + 4);
+ for (int i = 0; i < K; ++i)
+ {
+ bx(i) = point_set[i](0);
+ by(i) = point_set[i](1);
+ bz(i) = point_set[i](2);
+ }
+
+ for (int i = 0; i < 4; ++i)
+ {
+ bx(K + i) = start_end_derivative[i](0);
+ by(K + i) = start_end_derivative[i](1);
+ bz(K + i) = start_end_derivative[i](2);
+ }
+
+ // solve Ax = b
+ Eigen::VectorXd px = A.colPivHouseholderQr().solve(bx);
+ Eigen::VectorXd py = A.colPivHouseholderQr().solve(by);
+ Eigen::VectorXd pz = A.colPivHouseholderQr().solve(bz);
+
+ // convert to control pts
+ ctrl_pts.resize(3, K + 2);
+ ctrl_pts.row(0) = px.transpose();
+ ctrl_pts.row(1) = py.transpose();
+ ctrl_pts.row(2) = pz.transpose();
+
+ // cout << "[B-spline]: parameterization ok." << endl;
+ }
+
+ double UniformBspline::getTimeSum()
+ {
+ double tm, tmp;
+ if (getTimeSpan(tm, tmp))
+ return tmp - tm;
+ else
+ return -1.0;
+ }
+
+ double UniformBspline::getLength(const double &res)
+ {
+ double length = 0.0;
+ double dur = getTimeSum();
+ Eigen::VectorXd p_l = evaluateDeBoorT(0.0), p_n;
+ for (double t = res; t <= dur + 1e-4; t += res)
+ {
+ p_n = evaluateDeBoorT(t);
+ length += (p_n - p_l).norm();
+ p_l = p_n;
+ }
+ return length;
+ }
+
+ double UniformBspline::getJerk()
+ {
+ UniformBspline jerk_traj = getDerivative().getDerivative().getDerivative();
+
+ Eigen::VectorXd times = jerk_traj.getKnot();
+ Eigen::MatrixXd ctrl_pts = jerk_traj.getControlPoint();
+ int dimension = ctrl_pts.rows();
+
+ double jerk = 0.0;
+ for (int i = 0; i < ctrl_pts.cols(); ++i)
+ {
+ for (int j = 0; j < dimension; ++j)
+ {
+ jerk += (times(i + 1) - times(i)) * ctrl_pts(j, i) * ctrl_pts(j, i);
+ }
+ }
+
+ return jerk;
+ }
+
+ void UniformBspline::getMeanAndMaxVel(double &mean_v, double &max_v)
+ {
+ UniformBspline vel = getDerivative();
+ double tm, tmp;
+ vel.getTimeSpan(tm, tmp);
+
+ double max_vel = -1.0, mean_vel = 0.0;
+ int num = 0;
+ for (double t = tm; t <= tmp; t += 0.01)
+ {
+ Eigen::VectorXd vxd = vel.evaluateDeBoor(t);
+ double vn = vxd.norm();
+
+ mean_vel += vn;
+ ++num;
+ if (vn > max_vel)
+ {
+ max_vel = vn;
+ }
+ }
+
+ mean_vel = mean_vel / double(num);
+ mean_v = mean_vel;
+ max_v = max_vel;
+ }
+
+ void UniformBspline::getMeanAndMaxAcc(double &mean_a, double &max_a)
+ {
+ UniformBspline acc = getDerivative().getDerivative();
+ double tm, tmp;
+ acc.getTimeSpan(tm, tmp);
+
+ double max_acc = -1.0, mean_acc = 0.0;
+ int num = 0;
+ for (double t = tm; t <= tmp; t += 0.01)
+ {
+ Eigen::VectorXd axd = acc.evaluateDeBoor(t);
+ double an = axd.norm();
+
+ mean_acc += an;
+ ++num;
+ if (an > max_acc)
+ {
+ max_acc = an;
+ }
+ }
+
+ mean_acc = mean_acc / double(num);
+ mean_a = mean_acc;
+ max_a = max_acc;
+ }
+} // namespace ego_planner
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/CMakeLists.txt b/motion_planning/3d/ego_planner/planner/drone_detect/CMakeLists.txt
new file mode 100644
index 0000000..1eacae4
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/CMakeLists.txt
@@ -0,0 +1,130 @@
+cmake_minimum_required(VERSION 2.8.3)
+project(drone_detect)
+
+set(CMAKE_BUILD_TYPE "Release")
+ADD_COMPILE_OPTIONS(-std=c++11 )
+ADD_COMPILE_OPTIONS(-std=c++14 )
+## By adding -Wall and -Werror, the compiler does not ignore warnings anymore,
+## enforcing cleaner code.
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -Wall -g")
+
+## Find catkin macros and libraries
+find_package(catkin REQUIRED
+ COMPONENTS
+ roscpp
+ sensor_msgs
+ roslint
+ cv_bridge
+ message_filters
+)
+
+## Find system libraries
+find_package(Eigen3 REQUIRED)
+find_package(OpenCV REQUIRED)
+find_package(Boost REQUIRED)
+
+###################################
+## catkin specific configuration ##
+###################################
+## The catkin_package macro generates cmake config files for your package
+## Declare things to be passed to dependent projects
+## INCLUDE_DIRS: uncomment this if your package contains header files
+## LIBRARIES: libraries you create in this project that dependent projects also need
+## CATKIN_DEPENDS: catkin_packages dependent projects also need
+## DEPENDS: system dependencies of this project that dependent projects also need
+catkin_package(
+ INCLUDE_DIRS
+ include
+ ## This is only necessary because Eigen3 sets a non-standard EIGEN3_INCLUDE_DIR variable
+ ${EIGEN3_INCLUDE_DIR}
+ LIBRARIES
+ CATKIN_DEPENDS
+ roscpp
+ sensor_msgs
+ DEPENDS OpenCV Eigen Boost
+ ## find_package(Eigen3) provides a non standard EIGEN3_INCLUDE_DIR instead of Eigen3_INCLUDE_DIRS.
+ ## Therefore, the DEPEND does not work as expected and we need to add the directory to the INCLUDE_DIRS
+ # Eigen3
+
+ ## Boost is not part of the DEPENDS since it is only used in source files,
+ ## Dependees do not depend on Boost when they depend on this package.
+)
+
+###########
+## Build ##
+###########
+
+## Specify additional locations of header files
+## Your package locations should be listed before other locations
+include_directories(
+ include
+ ${catkin_INCLUDE_DIRS}
+ # Set manually because Eigen sets a non standard INCLUDE DIR
+ ${EIGEN3_INCLUDE_DIR}
+ # Set because Boost is an internal dependency, not transitive.
+ ${Boost_INCLUDE_DIRS}
+ ${OpenCV_INCLUDE_DIRS}
+
+)
+
+## Declare cpp executables
+add_executable(${PROJECT_NAME}
+ src/${PROJECT_NAME}_node.cpp
+ src/drone_detector.cpp
+)
+target_compile_features(${PROJECT_NAME} INTERFACE cxx_std_11)
+
+## Add dependencies to exported targets, like ROS msgs or srvs
+add_dependencies(${PROJECT_NAME}
+ ${catkin_EXPORTED_TARGETS}
+)
+
+
+target_link_libraries(${PROJECT_NAME}
+ ${catkin_LIBRARIES}
+ ${OpenCV_LIBS}
+)
+
+#############
+## Install ##
+#############
+
+# Mark executables and/or libraries for installation
+install(
+ TARGETS ${PROJECT_NAME}
+ ARCHIVE DESTINATION ${CATKIN_PACKAGE_LIB_DESTINATION}
+ LIBRARY DESTINATION ${CATKIN_PACKAGE_LIB_DESTINATION}
+ RUNTIME DESTINATION ${CATKIN_PACKAGE_BIN_DESTINATION}
+)
+
+# Mark cpp header files for installation
+install(
+ DIRECTORY include/${PROJECT_NAME}/
+ DESTINATION ${CATKIN_PACKAGE_INCLUDE_DESTINATION}
+ FILES_MATCHING PATTERN "*.h"
+)
+
+# Mark other files for installation
+install(
+ DIRECTORY doc
+ DESTINATION ${CATKIN_PACKAGE_SHARE_DESTINATION}
+)
+
+#############
+## Testing ##
+#############
+
+# if(${CATKIN_ENABLE_TESTING})
+# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pthread")
+# ## Add gtest based cpp test target and link libraries
+# catkin_add_gtest(${PROJECT_NAME}-test
+# test/test_drone_detector.cpp)
+
+# target_link_libraries(${PROJECT_NAME}-test)
+# endif()
+
+##########################
+## Static code analysis ##
+##########################
+
+roslint_cpp()
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/LICENSE b/motion_planning/3d/ego_planner/planner/drone_detect/LICENSE
new file mode 100644
index 0000000..9b02d46
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/LICENSE
@@ -0,0 +1,29 @@
+Software License Agreement (BSD License)
+
+Copyright (c) 2017, Peter Fankhauser
+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.
+
+* Neither the name of Autonomous Systems Lab nor ETH Zurich nor the names
+ of its contributors may be used to endorse or promote products derived from
+ this software without specific prior written permission.
+
+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/motion_planning/3d/ego_planner/planner/drone_detect/README.md b/motion_planning/3d/ego_planner/planner/drone_detect/README.md
new file mode 100644
index 0000000..2388e52
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/README.md
@@ -0,0 +1,97 @@
+# Drone Detect
+
+## Overview
+
+This is a package for detecting other drones in depth image and then calculating their true pose error in the world frame of the observer drone.
+
+
+
+
+
+
+
+## Usage
+
+You can launch the node alongside the main after you assigning the right topic name
+
+```
+roslaunch drone_detect drone_detect.launch
+```
+
+
+
+or add the following code in `run_in_sim.launch`
+
+```xml
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+```
+
+
+
+## Config files
+
+* **camera.yaml** : The camera intrinsics are stored in this file
+* **default.yaml**: Some parameters related to drone detection.
+
+```yaml
+debug_flag: true
+# the proportion of the pixel threshold to the total pixels projected from the drone to the depth map
+estimate/pixel_ratio: 0.1
+# the radius of the pose errror sphere
+estimate/max_pose_error: 0.4
+# the width and height of the model of drone
+estimate/drone_width: 0.5
+estimate/drone_height: 0.1
+```
+
+
+
+## Nodes
+
+
+#### Subscribed Topics
+
+* **`/depth`** ([sensor_msgs::Image])
+
+ The depth image from pcl_render_node.
+
+* **`/colordepth`**([sensor_msgs::Image])
+
+ The color image from pcl_render_node
+
+* **`/camera_pose`** ([geometry_msgs::PoseStamped])
+
+ The camere pose from pcl_render_node
+
+The above three topics are synchronized when in use, the callback function is **`rcvDepthColorCamPoseCallback`**
+
+- **`/dronex_odom_sub_`** ([nav_msgs::Odometry])
+
+ The odometry of other drones
+
+#### Published Topics
+
+* **`/new_depth`** ([sensor_msgs::Image])
+
+ The new depth image after erasing the moving drones
+
+* **`/new_colordepth`**([sensor_msgs::Image])
+
+ The color image with some debug marks
+
+* **`/camera_pose_error`** ([geometry_msgs::PoseStamped])
+
+ The pose error of detected drones in world frame of the observer drone.
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/config/camera.yaml b/motion_planning/3d/ego_planner/planner/drone_detect/config/camera.yaml
new file mode 100755
index 0000000..556e596
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/config/camera.yaml
@@ -0,0 +1,7 @@
+cam_width: 640
+cam_height: 480
+cam_fx: 386.02557373046875
+cam_fy: 386.02557373046875
+cam_cx: 321.8554382324219
+cam_cy: 241.2396240234375
+
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/config/default.yaml b/motion_planning/3d/ego_planner/planner/drone_detect/config/default.yaml
new file mode 100644
index 0000000..37a666d
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/config/default.yaml
@@ -0,0 +1,5 @@
+# debug_flag: true
+pixel_ratio: 0.1
+estimate/max_pose_error: 0.4
+estimate/drone_width: 0.5
+estimate/drone_height: 0.2
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/doc/demo.jpg b/motion_planning/3d/ego_planner/planner/drone_detect/doc/demo.jpg
new file mode 100644
index 0000000..690972d
Binary files /dev/null and b/motion_planning/3d/ego_planner/planner/drone_detect/doc/demo.jpg differ
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/doc/example.jpg b/motion_planning/3d/ego_planner/planner/drone_detect/doc/example.jpg
new file mode 100644
index 0000000..d047c4f
Binary files /dev/null and b/motion_planning/3d/ego_planner/planner/drone_detect/doc/example.jpg differ
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/include/drone_detector/drone_detector.h b/motion_planning/3d/ego_planner/planner/drone_detect/include/drone_detector/drone_detector.h
new file mode 100644
index 0000000..271210e
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/include/drone_detector/drone_detector.h
@@ -0,0 +1,156 @@
+#pragma once
+#include
+#include
+// ROS
+#include
+#include "std_msgs/String.h"
+#include "std_msgs/Bool.h"
+#include
+#include
+#include
+// synchronize topic
+#include
+#include
+#include
+#include
+
+#include
+
+//include opencv and eigen
+#include
+#include
+#include
+#include
+#include
+
+namespace detect {
+
+const int max_drone_num_ = 3;
+
+/*!
+ * Main class for the node to handle the ROS interfacing.
+ */
+class DroneDetector
+{
+ public:
+ /*!
+ * Constructor.
+ * @param nodeHandle the ROS node handle.
+ */
+ DroneDetector(ros::NodeHandle& nodeHandle);
+
+ /*!
+ * Destructor.
+ */
+ virtual ~DroneDetector();
+
+ void test();
+ private:
+ void readParameters();
+
+ // inline functions
+ double getDist2(const Eigen::Vector3d &p1, const Eigen::Vector3d &p2);
+ double getDist2(const Eigen::Vector4d &p1, const Eigen::Vector4d &p2);
+ Eigen::Vector4d depth2Pos(int u, int v, float depth);
+ Eigen::Vector4d depth2Pos(const Eigen::Vector2i &pixel, float depth) ;
+ Eigen::Vector2i pos2Depth(const Eigen::Vector4d &pose_in_camera) ;
+ bool isInSensorRange(const Eigen::Vector2i &pixel);
+
+ bool countPixel(int drone_id, Eigen::Vector2i &true_pixel, Eigen::Vector4d &true_pose_cam);
+ void detect(int drone_id, Eigen::Vector2i &true_pixel);
+
+ // subscribe callback function
+ void rcvDepthColorCamPoseCallback(const sensor_msgs::ImageConstPtr& depth_img, \
+ const sensor_msgs::ImageConstPtr& color_img, \
+ const geometry_msgs::PoseStampedConstPtr& camera_pose);
+
+ void rcvDepthCamPoseCallback(const sensor_msgs::ImageConstPtr& depth_img, \
+ const geometry_msgs::PoseStampedConstPtr& camera_pose);
+
+ void rcvMyOdomCallback(const nav_msgs::Odometry& odom);
+ void rcvDepthImgCallback(const sensor_msgs::ImageConstPtr& depth_img);
+
+
+ void rcvDroneOdomCallbackBase(const nav_msgs::Odometry& odom, const int drone_id);
+
+ void rcvDrone0OdomCallback(const nav_msgs::Odometry& odom);
+ void rcvDrone1OdomCallback(const nav_msgs::Odometry& odom);
+ void rcvDrone2OdomCallback(const nav_msgs::Odometry& odom);
+ void rcvDroneXOdomCallback(const nav_msgs::Odometry& odom);
+
+ //! ROS node handle.
+ ros::NodeHandle& nh_;
+
+ //! ROS topic subscriber.
+ // depth, colordepth, camera_pos subscriber
+ typedef message_filters::sync_policies::ApproximateTime SyncPolicyDepthColorImagePose;
+ typedef std::shared_ptr> SynchronizerDepthColorImagePose;
+ typedef message_filters::sync_policies::ApproximateTime SyncPolicyDepthImagePose;
+ typedef std::shared_ptr> SynchronizerDepthImagePose;
+
+ // std::shared_ptr> depth_img_sub_;
+ std::shared_ptr> colordepth_img_sub_;
+ std::shared_ptr> camera_pos_sub_;
+
+ SynchronizerDepthColorImagePose sync_depth_color_img_pose_;
+ SynchronizerDepthImagePose sync_depth_img_pose_;
+ // other drones subscriber
+ ros::Subscriber drone0_odom_sub_, drone1_odom_sub_, drone2_odom_sub_, droneX_odom_sub_;
+ std::string drone1_odom_topic_, drone2_odom_topic_;
+
+ ros::Subscriber my_odom_sub_, depth_img_sub_;
+ bool has_odom_;
+ nav_msgs::Odometry my_odom_;
+ // ROS topic publisher
+ // new_depth_img: erase the detected drones
+ // new_colordepth_img: for debug
+ ros::Publisher new_depth_img_pub_;
+ ros::Publisher debug_depth_img_pub_;
+
+ // parameters
+ //camera param
+ int img_width_, img_height_;
+ double fx_,fy_,cx_,cy_;
+
+ double max_pose_error_;
+ double max_pose_error2_;
+ double drone_width_, drone_height_;
+ double pixel_ratio_;
+ int pixel_threshold_;
+
+ // for debug
+ bool debug_flag_;
+ int debug_detect_result_[max_drone_num_];
+ std::stringstream debug_img_text_[max_drone_num_];
+ ros::Time debug_start_time_, debug_end_time_;
+
+ ros::Publisher debug_info_pub_;
+ ros::Publisher drone_pose_err_pub_[max_drone_num_];
+
+ int my_id_;
+ cv::Mat depth_img_, color_img_;
+
+ Eigen::Matrix4d cam2body_;
+ Eigen::Matrix4d cam2world_;
+ Eigen::Quaterniond cam2world_quat_;
+ Eigen::Vector4d my_pose_world_;
+ Eigen::Quaterniond my_attitude_world_;
+ Eigen::Vector4d my_last_pose_world_;
+ ros::Time my_last_odom_stamp_ = ros::TIME_MAX;
+ ros::Time my_last_camera_stamp_ = ros::TIME_MAX;
+
+ Eigen::Matrix4d drone2world_[max_drone_num_];
+ Eigen::Vector4d drone_pose_world_[max_drone_num_];
+ Eigen::Quaterniond drone_attitude_world_[max_drone_num_];
+ Eigen::Vector4d drone_pose_cam_[max_drone_num_];
+ Eigen::Vector2i drone_ref_pixel_[max_drone_num_];
+
+ std::vector hit_pixels_[max_drone_num_];
+ int valid_pixel_cnt_[max_drone_num_];
+
+ bool in_depth_[max_drone_num_] = {false};
+ cv::Point searchbox_lu_[max_drone_num_], searchbox_rd_[max_drone_num_];
+ cv::Point boundingbox_lu_[max_drone_num_], boundingbox_rd_[max_drone_num_];
+};
+
+} /* namespace */
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/launch/drone_detect.launch b/motion_planning/3d/ego_planner/planner/drone_detect/launch/drone_detect.launch
new file mode 100644
index 0000000..6aa4a25
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/launch/drone_detect.launch
@@ -0,0 +1,24 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/launch/ros_package_template_overlying_params.launch b/motion_planning/3d/ego_planner/planner/drone_detect/launch/ros_package_template_overlying_params.launch
new file mode 100644
index 0000000..0ac3e22
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/launch/ros_package_template_overlying_params.launch
@@ -0,0 +1,27 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/package.xml b/motion_planning/3d/ego_planner/planner/drone_detect/package.xml
new file mode 100644
index 0000000..bef3270
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/package.xml
@@ -0,0 +1,20 @@
+
+
+ drone_detect
+ 0.1.0
+ Detect other drones in depth image.
+ Jiangchao Zhu
+ BSD
+ Jiangchao Zhu
+
+
+ catkin
+
+ boost
+
+ eigen
+ roscpp
+ sensor_msgs
+
+ roslint
+
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/src/drone_detect_node.cpp b/motion_planning/3d/ego_planner/planner/drone_detect/src/drone_detect_node.cpp
new file mode 100644
index 0000000..8e29362
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/src/drone_detect_node.cpp
@@ -0,0 +1,14 @@
+#include
+#include "drone_detector/drone_detector.h"
+
+int main(int argc, char** argv)
+{
+ ros::init(argc, argv, "drone_detect");
+ ros::NodeHandle nh("~");
+
+ detect::DroneDetector drone_detector(nh);
+ drone_detector.test();
+
+ ros::spin();
+ return 0;
+}
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/src/drone_detector.cpp b/motion_planning/3d/ego_planner/planner/drone_detect/src/drone_detector.cpp
new file mode 100644
index 0000000..d3bacd3
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/src/drone_detector.cpp
@@ -0,0 +1,423 @@
+#include "drone_detector/drone_detector.h"
+
+// STD
+#include
+
+namespace detect {
+
+DroneDetector::DroneDetector(ros::NodeHandle& nodeHandle)
+ : nh_(nodeHandle)
+{
+ readParameters();
+
+ // depth_img_sub_.reset(new message_filters::Subscriber(nh_, "depth", 50, ros::TransportHints().tcpNoDelay()));
+ // colordepth_img_sub_.reset(new message_filters::Subscriber(nh_, "colordepth", 50));
+ // camera_pos_sub_.reset(new message_filters::Subscriber(nh_, "camera_pose", 50));
+
+ my_odom_sub_ = nh_.subscribe("odometry", 100, &DroneDetector::rcvMyOdomCallback, this, ros::TransportHints().tcpNoDelay());
+ depth_img_sub_ = nh_.subscribe("depth", 50, &DroneDetector::rcvDepthImgCallback, this, ros::TransportHints().tcpNoDelay());
+ // sync_depth_color_img_pose_->registerCallback(boost::bind(&DroneDetector::rcvDepthColorCamPoseCallback, this, _1, _2, _3));
+
+ // drone0_odom_sub_ = nh_.subscribe("drone0", 50, &DroneDetector::rcvDrone0OdomCallback, this);
+ // drone1_odom_sub_ = nh_.subscribe("drone1", 50, &DroneDetector::rcvDrone1OdomCallback, this);
+ // drone2_odom_sub_ = nh_.subscribe("drone2", 50, &DroneDetector::rcvDrone2OdomCallback, this);
+ droneX_odom_sub_ = nh_.subscribe("/others_odom", 100, &DroneDetector::rcvDroneXOdomCallback, this, ros::TransportHints().tcpNoDelay());
+
+ new_depth_img_pub_ = nh_.advertise("new_depth_image", 50);
+ debug_depth_img_pub_ = nh_.advertise("debug_depth_image", 50);
+
+ debug_info_pub_ = nh_.advertise("/debug_info", 50);
+
+ cam2body_ << 0.0, 0.0, 1.0, 0.0,
+ -1.0, 0.0, 0.0, 0.0,
+ 0.0, -1.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 1.0;
+
+ // init drone_pose_err_pub
+ for(int i = 0; i < max_drone_num_; i++) {
+ if(i != my_id_)
+ drone_pose_err_pub_[i] = nh_.advertise("drone"+std::to_string(i)+"to"+std::to_string(my_id_)+"_pose_err", 50);
+ }
+
+ ROS_INFO("Successfully launched node.");
+}
+
+DroneDetector::~DroneDetector()
+{
+}
+
+void DroneDetector::readParameters()
+{
+ // camera params
+ nh_.getParam("cam_width", img_width_);
+ nh_.getParam("cam_height", img_height_);
+ nh_.getParam("cam_fx", fx_);
+ nh_.getParam("cam_fy", fy_);
+ nh_.getParam("cam_cx", cx_);
+ nh_.getParam("cam_cy", cy_);
+
+ //
+ nh_.getParam("debug_flag", debug_flag_);
+ nh_.getParam("pixel_ratio", pixel_ratio_);
+ nh_.getParam("my_id", my_id_);
+ nh_.getParam("estimate/drone_width", drone_width_);
+ nh_.getParam("estimate/drone_height", drone_height_);
+ nh_.getParam("estimate/max_pose_error", max_pose_error_);
+
+ max_pose_error2_ = max_pose_error_*max_pose_error_;
+}
+
+// inline functions
+inline double DroneDetector::getDist2(const Eigen::Vector3d &p1, const Eigen::Vector3d &p2)
+{
+ double delta_x = p1(0)-p2(0);
+ double delta_y = p1(1)-p2(1);
+ double delta_z = p1(2)-p2(2);
+ return delta_x*delta_x+delta_y*delta_y+delta_z*delta_z;
+}
+
+inline double DroneDetector::getDist2(const Eigen::Vector4d &p1, const Eigen::Vector4d &p2)
+{
+ double delta_x = p1(0)-p2(0);
+ double delta_y = p1(1)-p2(1);
+ double delta_z = p1(2)-p2(2);
+ return delta_x*delta_x+delta_y*delta_y+delta_z*delta_z;
+}
+
+inline Eigen::Vector4d DroneDetector::depth2Pos(int u, int v, float depth)
+{
+ Eigen::Vector4d pose_in_camera;
+ pose_in_camera(0) = (u - cx_) * depth / fx_;
+ pose_in_camera(1) = (v - cy_) * depth / fy_;
+ pose_in_camera(2) = depth;
+ pose_in_camera(3) = 1.0;
+ return pose_in_camera;
+}
+
+inline Eigen::Vector4d DroneDetector::depth2Pos(const Eigen::Vector2i &pixel, float depth)
+{
+ Eigen::Vector4d pose_in_camera;
+ pose_in_camera(0) = (pixel(0) - cx_) * depth / fx_;
+ pose_in_camera(1) = (pixel(1) - cy_) * depth / fy_;
+ pose_in_camera(2) = depth;
+ pose_in_camera(3) = 1.0;
+ return pose_in_camera;
+}
+
+inline Eigen::Vector2i DroneDetector::pos2Depth(const Eigen::Vector4d &pose_in_camera)
+{
+ float depth = pose_in_camera(2);
+ Eigen::Vector2i pixel;
+ pixel(0) = pose_in_camera(0) * fx_ / depth + cx_ + 0.5;
+ pixel(1) = pose_in_camera(1) * fy_ / depth + cy_ + 0.5;
+ return pixel;
+}
+
+inline bool DroneDetector::isInSensorRange(const Eigen::Vector2i &pixel)
+{
+ if (pixel(0)>=0 && pixel(1) >= 0 && pixel(0) <= img_width_ && pixel(1) <= img_height_) return true;
+ else
+ return false;
+}
+
+void DroneDetector::rcvMyOdomCallback(const nav_msgs::Odometry& odom)
+{
+ my_odom_ = odom;
+ Eigen::Matrix4d body2world = Eigen::Matrix4d::Identity();
+
+ my_pose_world_(0) = odom.pose.pose.position.x;
+ my_pose_world_(1) = odom.pose.pose.position.y;
+ my_pose_world_(2) = odom.pose.pose.position.z;
+ my_pose_world_(3) = 1.0;
+ my_attitude_world_.x() = odom.pose.pose.orientation.x;
+ my_attitude_world_.y() = odom.pose.pose.orientation.y;
+ my_attitude_world_.z() = odom.pose.pose.orientation.z;
+ my_attitude_world_.w() = odom.pose.pose.orientation.w;
+ body2world.block<3,3>(0,0) = my_attitude_world_.toRotationMatrix();
+ body2world(0,3) = my_pose_world_(0);
+ body2world(1,3) = my_pose_world_(1);
+ body2world(2,3) = my_pose_world_(2);
+
+ //convert to cam pose
+ cam2world_ = body2world * cam2body_;
+ cam2world_quat_ = cam2world_.block<3,3>(0,0);
+
+ // my_last_odom_stamp_ = odom.header.stamp;
+
+ // my_last_pose_world_(0) = odom.pose.pose.position.x;
+ // my_last_pose_world_(1) = odom.pose.pose.position.y;
+ // my_last_pose_world_(2) = odom.pose.pose.position.z;
+ // my_last_pose_world_(3) = 1.0;
+
+ //publish tf
+ // static tf::TransformBroadcaster br;
+ // tf::Transform transform;
+ // transform.setOrigin( tf::Vector3(cam2world(0,3), cam2world(1,3), cam2world(2,3) ));
+ // transform.setRotation(tf::Quaternion(cam2world_quat.x(), cam2world_quat.y(), cam2world_quat.z(), cam2world_quat.w()));
+ // br.sendTransform(tf::StampedTransform(transform, my_last_odom_stamp, "world", "camera"));
+ //publish transform from world frame to quadrotor frame.
+}
+void DroneDetector::rcvDepthImgCallback(const sensor_msgs::ImageConstPtr& depth_img)
+{
+ /* get depth image */
+ cv_bridge::CvImagePtr cv_ptr;
+ cv_ptr = cv_bridge::toCvCopy(depth_img, depth_img->encoding);
+ cv_ptr->image.copyTo(depth_img_);
+
+ debug_start_time_ = ros::Time::now();
+
+ Eigen::Vector2i true_pixel[max_drone_num_];
+ for (int i = 0; i < max_drone_num_; i++) {
+ if (in_depth_[i]) {
+ detect(i, true_pixel[i]);
+ }
+ }
+
+ cv_bridge::CvImage out_msg;
+ for (int i = 0; i < max_drone_num_; i++) {
+ if (in_depth_[i]) {
+ // erase hit pixels in depth
+ for(int k = 0; k < int(hit_pixels_[i].size()); k++) {
+ // depth_img_.at(hit_pixels_[i][k](1), hit_pixels_[i][k](0)) = 0;
+ uint16_t *row_ptr;
+ row_ptr = depth_img_.ptr(hit_pixels_[i][k](1));
+ (*(row_ptr+hit_pixels_[i][k](0))) = 0.0;
+ }
+ }
+ }
+ debug_end_time_ = ros::Time::now();
+ // ROS_WARN("cost_total_time = %lf", (debug_end_time_ - debug_start_time_).toSec()*1000.0);
+ out_msg.header = depth_img->header;
+ out_msg.encoding = depth_img->encoding;
+ out_msg.image = depth_img_.clone();
+ new_depth_img_pub_.publish(out_msg.toImageMsg());
+
+ std_msgs::String msg;
+ std::stringstream ss;
+ if(debug_flag_) {
+ for (int i = 0; i < max_drone_num_; i++) {
+ if (in_depth_[i]) {
+ // add bound box in colormap
+ // cv::Rect rect(_bbox_lu.x, _bbox_lu.y, _bbox_rd.x, _bbox_rd.y);//左上坐标(x,y)和矩形的长(x)宽(y)
+ cv::rectangle(depth_img_, cv::Rect(searchbox_lu_[i], searchbox_rd_[i]), cv::Scalar(0, 0, 0), 5, cv::LINE_8, 0);
+ cv::rectangle(depth_img_, cv::Rect(boundingbox_lu_[i], boundingbox_rd_[i]), cv::Scalar(0, 0, 0), 5, cv::LINE_8, 0);
+ if (debug_detect_result_[i] == 1) {
+ ss << "no enough " << hit_pixels_[i].size();
+ } else if(debug_detect_result_[i] == 2) {
+ ss << "success";
+ }
+ } else {
+ ss << "no detect";
+ }
+ }
+ out_msg.header = depth_img->header;
+ out_msg.encoding = depth_img->encoding;
+ out_msg.image = depth_img_.clone();
+ debug_depth_img_pub_.publish(out_msg.toImageMsg());
+ msg.data = ss.str();
+ debug_info_pub_.publish(msg);
+ }
+}
+
+void DroneDetector::rcvDroneOdomCallbackBase(const nav_msgs::Odometry& odom, int drone_id)
+{
+ if (drone_id == my_id_) {
+ return;
+ }
+ Eigen::Matrix4d drone2world = Eigen::Matrix4d::Identity();
+ drone_pose_world_[drone_id](0) = odom.pose.pose.position.x;
+ drone_pose_world_[drone_id](1) = odom.pose.pose.position.y;
+ drone_pose_world_[drone_id](2) = odom.pose.pose.position.z;
+ drone_pose_world_[drone_id](3) = 1.0;
+
+ drone_attitude_world_[drone_id].x() = odom.pose.pose.orientation.x;
+ drone_attitude_world_[drone_id].y() = odom.pose.pose.orientation.y;
+ drone_attitude_world_[drone_id].z() = odom.pose.pose.orientation.z;
+ drone_attitude_world_[drone_id].w() = odom.pose.pose.orientation.w;
+ drone2world.block<3,3>(0,0) = drone_attitude_world_[drone_id].toRotationMatrix();
+
+ drone2world(0,3) = drone_pose_world_[drone_id](0);
+ drone2world(1,3) = drone_pose_world_[drone_id](1);
+ drone2world(2,3) = drone_pose_world_[drone_id](2);
+
+ drone_pose_cam_[drone_id] = cam2world_.inverse() * drone_pose_world_[drone_id];
+ // if the drone is in sensor range
+ drone_ref_pixel_[drone_id] = pos2Depth(drone_pose_cam_[drone_id]);
+ if (drone_pose_cam_[drone_id](2) > 0 && isInSensorRange(drone_ref_pixel_[drone_id])) {
+ in_depth_[drone_id] = true;
+ } else {
+ in_depth_[drone_id] = false;
+ debug_detect_result_[drone_id] = 0;
+ }
+}
+
+void DroneDetector::rcvDrone0OdomCallback(const nav_msgs::Odometry& odom)
+{
+ rcvDroneOdomCallbackBase(odom, 0);
+}
+
+void DroneDetector::rcvDrone1OdomCallback(const nav_msgs::Odometry& odom)
+{
+ rcvDroneOdomCallbackBase(odom, 1);
+}
+
+void DroneDetector::rcvDrone2OdomCallback(const nav_msgs::Odometry& odom)
+{
+ rcvDroneOdomCallbackBase(odom, 2);
+}
+
+void DroneDetector::rcvDroneXOdomCallback(const nav_msgs::Odometry& odom)
+{
+ std::string numstr = odom.child_frame_id.substr(6);
+ try
+ {
+ int drone_id = std::stoi(numstr);
+ rcvDroneOdomCallbackBase(odom, drone_id);
+ }
+ catch(const std::exception& e)
+ {
+ std::cout << e.what() << '\n';
+ }
+}
+
+bool DroneDetector::countPixel(int drone_id, Eigen::Vector2i &true_pixel, Eigen::Vector4d &true_pose_cam)
+{
+ boundingbox_lu_[drone_id].x = img_width_;
+ boundingbox_rd_[drone_id].x = 0;
+ boundingbox_lu_[drone_id].y = img_height_;
+ boundingbox_rd_[drone_id].y = 0;
+
+ valid_pixel_cnt_[drone_id] = 0;
+ hit_pixels_[drone_id].clear();
+
+ Eigen::Vector2i tmp_pixel;
+ Eigen::Vector4d tmp_pose_cam;
+ int search_radius = 2*max_pose_error_*fx_/drone_pose_cam_[drone_id](2);
+ float depth;
+ searchbox_lu_[drone_id].x = drone_ref_pixel_[drone_id](0) - search_radius;
+ searchbox_lu_[drone_id].y = drone_ref_pixel_[drone_id](1) - search_radius;
+ searchbox_rd_[drone_id].x = drone_ref_pixel_[drone_id](0) + search_radius;
+ searchbox_rd_[drone_id].y = drone_ref_pixel_[drone_id](1) + search_radius;
+ // check the tmp_p around ref_pixel
+ for(int i = -search_radius; i <= search_radius; i++)
+ for(int j = -search_radius; j <= search_radius; j++)
+ {
+ tmp_pixel(0) = drone_ref_pixel_[drone_id](0) + j;
+ tmp_pixel(1) = drone_ref_pixel_[drone_id](1) + i;
+ if(tmp_pixel(0) < 0 || tmp_pixel(0) >= img_width_ || tmp_pixel(1) < 0 || tmp_pixel(1) >= img_height_)
+ continue;
+ // depth = depth_img_.at(tmp_pixel(1), tmp_pixel(0));
+ uint16_t *row_ptr;
+ row_ptr = depth_img_.ptr(tmp_pixel(1));
+ depth = (*(row_ptr+tmp_pixel(0))) / 1000.0;
+ // ROS_WARN("depth = %lf", depth);
+ // get tmp_pose in cam frame
+ tmp_pose_cam = depth2Pos(tmp_pixel(0), tmp_pixel(1), depth);
+ double dist2 = getDist2(tmp_pose_cam, drone_pose_cam_[drone_id]);
+ // ROS_WARN("dist2 = %lf", dist2);
+ if (dist2 < max_pose_error2_) {
+ valid_pixel_cnt_[drone_id]++;
+ hit_pixels_[drone_id].push_back(tmp_pixel);
+ boundingbox_lu_[drone_id].x = tmp_pixel(0) < boundingbox_lu_[drone_id].x ? tmp_pixel(0) : boundingbox_lu_[drone_id].x;
+ boundingbox_lu_[drone_id].y = tmp_pixel(1) < boundingbox_lu_[drone_id].y ? tmp_pixel(1) : boundingbox_lu_[drone_id].y;
+ boundingbox_rd_[drone_id].x = tmp_pixel(0) > boundingbox_rd_[drone_id].x ? tmp_pixel(0) : boundingbox_rd_[drone_id].x;
+ boundingbox_rd_[drone_id].y = tmp_pixel(1) > boundingbox_rd_[drone_id].y ? tmp_pixel(1) : boundingbox_rd_[drone_id].y;
+ }
+ }
+ pixel_threshold_ = (drone_width_*fx_/drone_pose_cam_[drone_id](2)) * (drone_height_*fy_/drone_pose_cam_[drone_id](2))*pixel_ratio_;
+ if (valid_pixel_cnt_[drone_id] > pixel_threshold_) {
+ int step = 1, size = (boundingbox_rd_[drone_id].y-boundingbox_lu_[drone_id].y) < (boundingbox_rd_[drone_id].x-boundingbox_lu_[drone_id].x) ? (boundingbox_rd_[drone_id].y-boundingbox_lu_[drone_id].y) : (boundingbox_rd_[drone_id].x-boundingbox_lu_[drone_id].x);
+ int init_x = (boundingbox_lu_[drone_id].x+boundingbox_rd_[drone_id].x)/2, init_y = (boundingbox_lu_[drone_id].y+boundingbox_rd_[drone_id].y)/2;
+ int x_flag = 1, y_flag = 1;
+ int x_idx = 0, y_idx = 0;
+ uint16_t *row_ptr;
+ row_ptr = depth_img_.ptr(tmp_pixel(1));
+ depth = (*(row_ptr+tmp_pixel(0))) / 1000.0;
+ tmp_pose_cam = depth2Pos(init_x, init_y, depth);
+ if (getDist2(tmp_pose_cam, drone_pose_cam_[drone_id]) < max_pose_error2_){
+ true_pixel(0) = init_x;
+ true_pixel(1) = init_y;
+ true_pose_cam = tmp_pose_cam;
+ return true;
+ }
+ while(step(tmp_pixel(1));
+ depth = (*(row_ptr+tmp_pixel(0))) / 1000.0;
+ tmp_pose_cam = depth2Pos(init_x, init_y, depth);
+ if (getDist2(tmp_pose_cam, drone_pose_cam_[drone_id]) < max_pose_error2_) {
+ true_pixel(0) = init_x;
+ true_pixel(1) = init_y;
+ true_pose_cam = tmp_pose_cam;
+ return true;
+ }
+ x_idx++;
+ }
+ x_idx = 0;
+ x_flag = -x_flag;
+ while(y_idx(tmp_pixel(1));
+ depth = (*(row_ptr+tmp_pixel(0))) / 1000.0;
+ tmp_pose_cam = depth2Pos(init_x, init_y, depth);
+ if (getDist2(tmp_pose_cam, drone_pose_cam_[drone_id]) < max_pose_error2_){
+ true_pixel(0) = init_x;
+ true_pixel(1) = init_y;
+ true_pose_cam = tmp_pose_cam;
+ return true;
+ }
+ y_idx++;
+ }
+ y_idx = 0;
+ y_flag = -y_flag;
+ step++;
+ }
+ while(x_idx(tmp_pixel(1));
+ depth = (*(row_ptr+tmp_pixel(0))) / 1000.0;
+ tmp_pose_cam = depth2Pos(init_x, init_y, depth);
+ if (getDist2(tmp_pose_cam, drone_pose_cam_[drone_id]) < max_pose_error2_){
+ true_pixel(0) = init_x;
+ true_pixel(1) = init_y;
+ true_pose_cam = tmp_pose_cam;
+ return true;
+ }
+ x_idx++;
+ }
+ }
+ return false;
+}
+
+void DroneDetector::detect(int drone_id, Eigen::Vector2i &true_pixel)
+{
+ Eigen::Vector4d true_pose_cam, pose_error;
+ bool found = countPixel(drone_id, true_pixel, true_pose_cam);
+ if (found) {
+ // ROS_WARN("FOUND");
+ pose_error = cam2world_*true_pose_cam - drone_pose_world_[drone_id];
+ debug_detect_result_[drone_id] = 2;
+
+ geometry_msgs::PoseStamped out_msg;
+ out_msg.header.stamp = my_last_camera_stamp_;
+ out_msg.header.frame_id = "/drone_detect";
+ out_msg.pose.position.x = pose_error(0);
+ out_msg.pose.position.y = pose_error(1);
+ out_msg.pose.position.z = pose_error(2);
+ drone_pose_err_pub_[drone_id].publish(out_msg);
+
+ } else {
+ // ROS_WARN("NOT FOUND");
+ debug_detect_result_[drone_id] = 1;
+ }
+}
+
+void DroneDetector::test() {
+ ROS_WARN("my_id = %d", my_id_);
+}
+
+} /* namespace */
diff --git a/motion_planning/3d/ego_planner/planner/drone_detect/test/test_drone_detector.cpp b/motion_planning/3d/ego_planner/planner/drone_detect/test/test_drone_detector.cpp
new file mode 100644
index 0000000..2fc9dd5
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/drone_detect/test/test_drone_detector.cpp
@@ -0,0 +1,10 @@
+// gtest
+#include
+
+// Run all the tests that were declared with TEST()
+int main(int argc, char **argv)
+{
+ testing::InitGoogleTest(&argc, argv);
+ srand((int)time(0));
+ return RUN_ALL_TESTS();
+}
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/path_searching/CMakeLists.txt b/motion_planning/3d/ego_planner/planner/path_searching/CMakeLists.txt
new file mode 100755
index 0000000..13f1f66
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/path_searching/CMakeLists.txt
@@ -0,0 +1,42 @@
+cmake_minimum_required(VERSION 2.8.3)
+project(path_searching)
+
+set(CMAKE_BUILD_TYPE "Release")
+ADD_COMPILE_OPTIONS(-std=c++11 )
+ADD_COMPILE_OPTIONS(-std=c++14 )
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -Wall -g")
+
+find_package(Eigen3 REQUIRED)
+find_package(PCL 1.7 REQUIRED)
+
+find_package(catkin REQUIRED COMPONENTS
+ roscpp
+ rospy
+ std_msgs
+ visualization_msgs
+ plan_env
+ cv_bridge
+)
+
+
+catkin_package(
+ INCLUDE_DIRS include
+ LIBRARIES path_searching
+ CATKIN_DEPENDS plan_env
+# DEPENDS system_lib
+)
+
+include_directories(
+ SYSTEM
+ include
+ ${catkin_INCLUDE_DIRS}
+ ${Eigen3_INCLUDE_DIRS}
+ ${PCL_INCLUDE_DIRS}
+)
+
+add_library( path_searching
+ src/dyn_a_star.cpp
+ )
+target_link_libraries( path_searching
+ ${catkin_LIBRARIES}
+ )
diff --git a/motion_planning/3d/ego_planner/planner/path_searching/include/path_searching/dyn_a_star.h b/motion_planning/3d/ego_planner/planner/path_searching/include/path_searching/dyn_a_star.h
new file mode 100755
index 0000000..388668b
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/path_searching/include/path_searching/dyn_a_star.h
@@ -0,0 +1,115 @@
+#ifndef _DYN_A_STAR_H_
+#define _DYN_A_STAR_H_
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+constexpr double inf = 1 >> 20;
+struct GridNode;
+typedef GridNode *GridNodePtr;
+
+struct GridNode
+{
+ enum enum_state
+ {
+ OPENSET = 1,
+ CLOSEDSET = 2,
+ UNDEFINED = 3
+ };
+
+ int rounds{0}; // Distinguish every call
+ enum enum_state state
+ {
+ UNDEFINED
+ };
+ Eigen::Vector3i index;
+
+ double gScore{inf}, fScore{inf};
+ GridNodePtr cameFrom{NULL};
+};
+
+class NodeComparator
+{
+public:
+ bool operator()(GridNodePtr node1, GridNodePtr node2)
+ {
+ return node1->fScore > node2->fScore;
+ }
+};
+
+class AStar
+{
+private:
+ GridMap::Ptr grid_map_;
+
+ inline void coord2gridIndexFast(const double x, const double y, const double z, int &id_x, int &id_y, int &id_z);
+
+ double getDiagHeu(GridNodePtr node1, GridNodePtr node2);
+ double getManhHeu(GridNodePtr node1, GridNodePtr node2);
+ double getEuclHeu(GridNodePtr node1, GridNodePtr node2);
+ inline double getHeu(GridNodePtr node1, GridNodePtr node2);
+
+ bool ConvertToIndexAndAdjustStartEndPoints(const Eigen::Vector3d start_pt, const Eigen::Vector3d end_pt, Eigen::Vector3i &start_idx, Eigen::Vector3i &end_idx);
+
+ inline Eigen::Vector3d Index2Coord(const Eigen::Vector3i &index) const;
+ inline bool Coord2Index(const Eigen::Vector3d &pt, Eigen::Vector3i &idx) const;
+
+ //bool (*checkOccupancyPtr)( const Eigen::Vector3d &pos );
+
+ inline bool checkOccupancy(const Eigen::Vector3d &pos) { return (bool)grid_map_->getInflateOccupancy(pos); }
+
+ std::vector retrievePath(GridNodePtr current);
+
+ double step_size_, inv_step_size_;
+ Eigen::Vector3d center_;
+ Eigen::Vector3i CENTER_IDX_, POOL_SIZE_;
+ const double tie_breaker_ = 1.0 + 1.0 / 10000;
+
+ std::vector gridPath_;
+
+ GridNodePtr ***GridNodeMap_;
+ std::priority_queue, NodeComparator> openSet_;
+
+ int rounds_{0};
+
+public:
+ typedef std::shared_ptr Ptr;
+
+ AStar(){};
+ ~AStar();
+
+ void initGridMap(GridMap::Ptr occ_map, const Eigen::Vector3i pool_size);
+
+ bool AstarSearch(const double step_size, Eigen::Vector3d start_pt, Eigen::Vector3d end_pt);
+
+ std::vector getPath();
+};
+
+inline double AStar::getHeu(GridNodePtr node1, GridNodePtr node2)
+{
+ return tie_breaker_ * getDiagHeu(node1, node2);
+}
+
+inline Eigen::Vector3d AStar::Index2Coord(const Eigen::Vector3i &index) const
+{
+ return ((index - CENTER_IDX_).cast() * step_size_) + center_;
+};
+
+inline bool AStar::Coord2Index(const Eigen::Vector3d &pt, Eigen::Vector3i &idx) const
+{
+ idx = ((pt - center_) * inv_step_size_ + Eigen::Vector3d(0.5, 0.5, 0.5)).cast() + CENTER_IDX_;
+
+ if (idx(0) < 0 || idx(0) >= POOL_SIZE_(0) || idx(1) < 0 || idx(1) >= POOL_SIZE_(1) || idx(2) < 0 || idx(2) >= POOL_SIZE_(2))
+ {
+ ROS_ERROR("Ran out of pool, index=%d %d %d", idx(0), idx(1), idx(2));
+ return false;
+ }
+
+ return true;
+};
+
+#endif
diff --git a/motion_planning/3d/ego_planner/planner/path_searching/package.xml b/motion_planning/3d/ego_planner/planner/path_searching/package.xml
new file mode 100755
index 0000000..0394974
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/path_searching/package.xml
@@ -0,0 +1,71 @@
+
+
+ path_searching
+ 0.0.0
+ The path_searching package
+
+
+
+
+ iszhouxin
+
+
+
+
+
+ TODO
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ catkin
+ roscpp
+ rospy
+ std_msgs
+ plan_env
+ roscpp
+ rospy
+ std_msgs
+ plan_env
+ roscpp
+ rospy
+ std_msgs
+ plan_env
+
+
+
+
+
+
+
+
diff --git a/motion_planning/3d/ego_planner/planner/path_searching/src/dyn_a_star.cpp b/motion_planning/3d/ego_planner/planner/path_searching/src/dyn_a_star.cpp
new file mode 100755
index 0000000..744daa9
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/path_searching/src/dyn_a_star.cpp
@@ -0,0 +1,260 @@
+#include "path_searching/dyn_a_star.h"
+
+using namespace std;
+using namespace Eigen;
+
+AStar::~AStar()
+{
+ for (int i = 0; i < POOL_SIZE_(0); i++)
+ for (int j = 0; j < POOL_SIZE_(1); j++)
+ for (int k = 0; k < POOL_SIZE_(2); k++)
+ delete GridNodeMap_[i][j][k];
+}
+
+void AStar::initGridMap(GridMap::Ptr occ_map, const Eigen::Vector3i pool_size)
+{
+ POOL_SIZE_ = pool_size;
+ CENTER_IDX_ = pool_size / 2;
+
+ GridNodeMap_ = new GridNodePtr **[POOL_SIZE_(0)];
+ for (int i = 0; i < POOL_SIZE_(0); i++)
+ {
+ GridNodeMap_[i] = new GridNodePtr *[POOL_SIZE_(1)];
+ for (int j = 0; j < POOL_SIZE_(1); j++)
+ {
+ GridNodeMap_[i][j] = new GridNodePtr[POOL_SIZE_(2)];
+ for (int k = 0; k < POOL_SIZE_(2); k++)
+ {
+ GridNodeMap_[i][j][k] = new GridNode;
+ }
+ }
+ }
+
+ grid_map_ = occ_map;
+}
+
+double AStar::getDiagHeu(GridNodePtr node1, GridNodePtr node2)
+{
+ double dx = abs(node1->index(0) - node2->index(0));
+ double dy = abs(node1->index(1) - node2->index(1));
+ double dz = abs(node1->index(2) - node2->index(2));
+
+ double h = 0.0;
+ int diag = min(min(dx, dy), dz);
+ dx -= diag;
+ dy -= diag;
+ dz -= diag;
+
+ if (dx == 0)
+ {
+ h = 1.0 * sqrt(3.0) * diag + sqrt(2.0) * min(dy, dz) + 1.0 * abs(dy - dz);
+ }
+ if (dy == 0)
+ {
+ h = 1.0 * sqrt(3.0) * diag + sqrt(2.0) * min(dx, dz) + 1.0 * abs(dx - dz);
+ }
+ if (dz == 0)
+ {
+ h = 1.0 * sqrt(3.0) * diag + sqrt(2.0) * min(dx, dy) + 1.0 * abs(dx - dy);
+ }
+ return h;
+}
+
+double AStar::getManhHeu(GridNodePtr node1, GridNodePtr node2)
+{
+ double dx = abs(node1->index(0) - node2->index(0));
+ double dy = abs(node1->index(1) - node2->index(1));
+ double dz = abs(node1->index(2) - node2->index(2));
+
+ return dx + dy + dz;
+}
+
+double AStar::getEuclHeu(GridNodePtr node1, GridNodePtr node2)
+{
+ return (node2->index - node1->index).norm();
+}
+
+vector AStar::retrievePath(GridNodePtr current)
+{
+ vector path;
+ path.push_back(current);
+
+ while (current->cameFrom != NULL)
+ {
+ current = current->cameFrom;
+ path.push_back(current);
+ }
+
+ return path;
+}
+
+bool AStar::ConvertToIndexAndAdjustStartEndPoints(Vector3d start_pt, Vector3d end_pt, Vector3i &start_idx, Vector3i &end_idx)
+{
+ if (!Coord2Index(start_pt, start_idx) || !Coord2Index(end_pt, end_idx))
+ return false;
+
+ if (checkOccupancy(Index2Coord(start_idx)))
+ {
+ //ROS_WARN("Start point is insdide an obstacle.");
+ do
+ {
+ start_pt = (start_pt - end_pt).normalized() * step_size_ + start_pt;
+ if (!Coord2Index(start_pt, start_idx))
+ return false;
+ } while (checkOccupancy(Index2Coord(start_idx)));
+ }
+
+ if (checkOccupancy(Index2Coord(end_idx)))
+ {
+ //ROS_WARN("End point is insdide an obstacle.");
+ do
+ {
+ end_pt = (end_pt - start_pt).normalized() * step_size_ + end_pt;
+ if (!Coord2Index(end_pt, end_idx))
+ return false;
+ } while (checkOccupancy(Index2Coord(end_idx)));
+ }
+
+ return true;
+}
+
+bool AStar::AstarSearch(const double step_size, Vector3d start_pt, Vector3d end_pt)
+{
+ ros::Time time_1 = ros::Time::now();
+ ++rounds_;
+
+ step_size_ = step_size;
+ inv_step_size_ = 1 / step_size;
+ center_ = (start_pt + end_pt) / 2;
+
+ Vector3i start_idx, end_idx;
+ if (!ConvertToIndexAndAdjustStartEndPoints(start_pt, end_pt, start_idx, end_idx))
+ {
+ ROS_ERROR("Unable to handle the initial or end point, force return!");
+ return false;
+ }
+
+ // if ( start_pt(0) > -1 && start_pt(0) < 0 )
+ // cout << "start_pt=" << start_pt.transpose() << " end_pt=" << end_pt.transpose() << endl;
+
+ GridNodePtr startPtr = GridNodeMap_[start_idx(0)][start_idx(1)][start_idx(2)];
+ GridNodePtr endPtr = GridNodeMap_[end_idx(0)][end_idx(1)][end_idx(2)];
+
+ std::priority_queue, NodeComparator> empty;
+ openSet_.swap(empty);
+
+ GridNodePtr neighborPtr = NULL;
+ GridNodePtr current = NULL;
+
+ startPtr->index = start_idx;
+ startPtr->rounds = rounds_;
+ startPtr->gScore = 0;
+ startPtr->fScore = getHeu(startPtr, endPtr);
+ startPtr->state = GridNode::OPENSET; //put start node in open set
+ startPtr->cameFrom = NULL;
+ openSet_.push(startPtr); //put start in open set
+
+ endPtr->index = end_idx;
+
+ double tentative_gScore;
+
+ int num_iter = 0;
+ while (!openSet_.empty())
+ {
+ num_iter++;
+ current = openSet_.top();
+ openSet_.pop();
+
+ // if ( num_iter < 10000 )
+ // cout << "current=" << current->index.transpose() << endl;
+
+ if (current->index(0) == endPtr->index(0) && current->index(1) == endPtr->index(1) && current->index(2) == endPtr->index(2))
+ {
+ // ros::Time time_2 = ros::Time::now();
+ // printf("\033[34mA star iter:%d, time:%.3f\033[0m\n",num_iter, (time_2 - time_1).toSec()*1000);
+ // if((time_2 - time_1).toSec() > 0.1)
+ // ROS_WARN("Time consume in A star path finding is %f", (time_2 - time_1).toSec() );
+ gridPath_ = retrievePath(current);
+ return true;
+ }
+ current->state = GridNode::CLOSEDSET; //move current node from open set to closed set.
+
+ for (int dx = -1; dx <= 1; dx++)
+ for (int dy = -1; dy <= 1; dy++)
+ for (int dz = -1; dz <= 1; dz++)
+ {
+ if (dx == 0 && dy == 0 && dz == 0)
+ continue;
+
+ Vector3i neighborIdx;
+ neighborIdx(0) = (current->index)(0) + dx;
+ neighborIdx(1) = (current->index)(1) + dy;
+ neighborIdx(2) = (current->index)(2) + dz;
+
+ if (neighborIdx(0) < 1 || neighborIdx(0) >= POOL_SIZE_(0) - 1 || neighborIdx(1) < 1 || neighborIdx(1) >= POOL_SIZE_(1) - 1 || neighborIdx(2) < 1 || neighborIdx(2) >= POOL_SIZE_(2) - 1)
+ {
+ continue;
+ }
+
+ neighborPtr = GridNodeMap_[neighborIdx(0)][neighborIdx(1)][neighborIdx(2)];
+ neighborPtr->index = neighborIdx;
+
+ bool flag_explored = neighborPtr->rounds == rounds_;
+
+ if (flag_explored && neighborPtr->state == GridNode::CLOSEDSET)
+ {
+ continue; //in closed set.
+ }
+
+ neighborPtr->rounds = rounds_;
+
+ if (checkOccupancy(Index2Coord(neighborPtr->index)))
+ {
+ continue;
+ }
+
+ double static_cost = sqrt(dx * dx + dy * dy + dz * dz);
+ tentative_gScore = current->gScore + static_cost;
+
+ if (!flag_explored)
+ {
+ //discover a new node
+ neighborPtr->state = GridNode::OPENSET;
+ neighborPtr->cameFrom = current;
+ neighborPtr->gScore = tentative_gScore;
+ neighborPtr->fScore = tentative_gScore + getHeu(neighborPtr, endPtr);
+ openSet_.push(neighborPtr); //put neighbor in open set and record it.
+ }
+ else if (tentative_gScore < neighborPtr->gScore)
+ { //in open set and need update
+ neighborPtr->cameFrom = current;
+ neighborPtr->gScore = tentative_gScore;
+ neighborPtr->fScore = tentative_gScore + getHeu(neighborPtr, endPtr);
+ }
+ }
+ ros::Time time_2 = ros::Time::now();
+ if ((time_2 - time_1).toSec() > 0.2)
+ {
+ ROS_WARN("Failed in A star path searching !!! 0.2 seconds time limit exceeded.");
+ return false;
+ }
+ }
+
+ ros::Time time_2 = ros::Time::now();
+
+ if ((time_2 - time_1).toSec() > 0.1)
+ ROS_WARN("Time consume in A star path finding is %.3fs, iter=%d", (time_2 - time_1).toSec(), num_iter);
+
+ return false;
+}
+
+vector AStar::getPath()
+{
+ vector path;
+
+ for (auto ptr : gridPath_)
+ path.push_back(Index2Coord(ptr->index));
+
+ reverse(path.begin(), path.end());
+ return path;
+}
diff --git a/motion_planning/3d/ego_planner/planner/plan_env/CMakeLists.txt b/motion_planning/3d/ego_planner/planner/plan_env/CMakeLists.txt
new file mode 100755
index 0000000..1567f41
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/plan_env/CMakeLists.txt
@@ -0,0 +1,57 @@
+cmake_minimum_required(VERSION 2.8.3)
+project(plan_env)
+
+set(CMAKE_BUILD_TYPE "Release")
+ADD_COMPILE_OPTIONS(-std=c++11 )
+ADD_COMPILE_OPTIONS(-std=c++14 )
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -Wall -g")
+
+find_package(OpenCV REQUIRED)
+
+
+find_package(catkin REQUIRED COMPONENTS
+ roscpp
+ rospy
+ std_msgs
+ visualization_msgs
+ cv_bridge
+ message_filters
+)
+
+find_package(Eigen3 REQUIRED)
+find_package(PCL 1.7 REQUIRED)
+
+catkin_package(
+ INCLUDE_DIRS include
+ LIBRARIES plan_env
+ CATKIN_DEPENDS roscpp std_msgs
+# DEPENDS system_lib
+)
+
+include_directories(
+ SYSTEM
+ include
+ ${catkin_INCLUDE_DIRS}
+ ${Eigen3_INCLUDE_DIRS}
+ ${PCL_INCLUDE_DIRS}
+ ${OpenCV_INCLUDE_DIRS}
+)
+
+link_directories(${PCL_LIBRARY_DIRS})
+
+add_library( plan_env
+ src/grid_map.cpp
+ src/raycast.cpp
+ src/obj_predictor.cpp
+ )
+target_link_libraries( plan_env
+ ${catkin_LIBRARIES}
+ ${PCL_LIBRARIES}
+ )
+
+add_executable(obj_generator
+ src/obj_generator.cpp
+)
+target_link_libraries(obj_generator
+ ${catkin_LIBRARIES}
+ )
diff --git a/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/grid_map.h b/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/grid_map.h
new file mode 100644
index 0000000..bba2603
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/grid_map.h
@@ -0,0 +1,804 @@
+#ifndef _GRID_MAP_H
+#define _GRID_MAP_H
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+#include
+
+#define logit(x) (log((x) / (1 - (x))))
+
+using namespace std;
+
+// voxel hashing
+template
+struct matrix_hash : std::unary_function {
+ std::size_t operator()(T const& matrix) const {
+ size_t seed = 0;
+ for (size_t i = 0; i < matrix.size(); ++i) {
+ auto elem = *(matrix.data() + i);
+ seed ^= std::hash()(elem) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+ }
+ return seed;
+ }
+};
+
+// constant parameters
+
+struct MappingParameters {
+
+ /* map properties */
+ Eigen::Vector3d map_origin_, map_size_;
+ Eigen::Vector3d map_min_boundary_, map_max_boundary_; // map range in pos
+ Eigen::Vector3i map_voxel_num_; // map range in index
+ Eigen::Vector3d local_update_range_;
+ double resolution_, resolution_inv_;
+ double obstacles_inflation_;
+ string frame_id_;
+ int pose_type_;
+
+ /* camera parameters */
+ double cx_, cy_, fx_, fy_;
+
+ /* time out */
+ double odom_depth_timeout_;
+
+ /* depth image projection filtering */
+ double depth_filter_maxdist_, depth_filter_mindist_, depth_filter_tolerance_;
+ int depth_filter_margin_;
+ bool use_depth_filter_;
+ double k_depth_scaling_factor_;
+ int skip_pixel_;
+
+ /* raycasting */
+ double p_hit_, p_miss_, p_min_, p_max_, p_occ_; // occupancy probability
+ double prob_hit_log_, prob_miss_log_, clamp_min_log_, clamp_max_log_,
+ min_occupancy_log_; // logit of occupancy probability
+ double min_ray_length_, max_ray_length_; // range of doing raycasting
+
+ /* local map update and clear */
+ int local_map_margin_;
+
+ /* visualization and computation time display */
+ double visualization_truncate_height_, virtual_ceil_height_, ground_height_, virtual_ceil_yp_, virtual_ceil_yn_;
+ bool show_occ_time_;
+
+ /* active mapping */
+ double unknown_flag_;
+};
+
+// intermediate mapping data for fusion
+
+struct MappingData {
+ // main map data, occupancy of each voxel and Euclidean distance
+
+ std::vector occupancy_buffer_;
+ std::vector occupancy_buffer_inflate_;
+
+ // camera position and pose data
+
+ Eigen::Vector3d camera_pos_, last_camera_pos_;
+ Eigen::Matrix3d camera_r_m_, last_camera_r_m_;
+ Eigen::Matrix4d cam2body_;
+
+ // depth image data
+
+ cv::Mat depth_image_, last_depth_image_;
+ int image_cnt_;
+
+ // flags of map state
+
+ bool occ_need_update_, local_updated_;
+ bool has_first_depth_;
+ bool has_odom_, has_cloud_;
+
+ // odom_depth_timeout_
+ ros::Time last_occ_update_time_;
+ bool flag_depth_odom_timeout_;
+ bool flag_use_depth_fusion;
+
+ // depth image projected point cloud
+
+ vector proj_points_;
+ int proj_points_cnt;
+
+ // flag buffers for speeding up raycasting
+
+ vector count_hit_, count_hit_and_miss_;
+ vector flag_traverse_, flag_rayend_;
+ char raycast_num_;
+ queue cache_voxel_;
+
+ // range of updating grid
+
+ Eigen::Vector3i local_bound_min_, local_bound_max_;
+
+ // computation time
+
+ double fuse_time_, max_fuse_time_;
+ int update_num_;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+};
+
+class GridMap {
+public:
+ GridMap() {}
+ ~GridMap() {}
+
+ enum { POSE_STAMPED = 1, ODOMETRY = 2, INVALID_IDX = -10000 };
+
+ // occupancy map management
+ void resetBuffer();
+ void resetBuffer(Eigen::Vector3d min, Eigen::Vector3d max);
+
+ inline void posToIndex(const Eigen::Vector3d& pos, Eigen::Vector3i& id);
+ inline void indexToPos(const Eigen::Vector3i& id, Eigen::Vector3d& pos);
+ inline int toAddress(const Eigen::Vector3i& id);
+ inline int toAddress(int& x, int& y, int& z);
+ inline bool isInMap(const Eigen::Vector3d& pos);
+ inline bool isInMap(const Eigen::Vector3i& idx);
+
+ inline void setOccupancy(Eigen::Vector3d pos, double occ = 1);
+ inline void setOccupied(Eigen::Vector3d pos);
+ inline int getOccupancy(Eigen::Vector3d pos);
+ inline int getOccupancy(Eigen::Vector3i id);
+ inline int getInflateOccupancy(Eigen::Vector3d pos);
+
+ inline void boundIndex(Eigen::Vector3i& id);
+ inline bool isUnknown(const Eigen::Vector3i& id);
+ inline bool isUnknown(const Eigen::Vector3d& pos);
+ inline bool isKnownFree(const Eigen::Vector3i& id);
+ inline bool isKnownOccupied(const Eigen::Vector3i& id);
+
+ void initMap(ros::NodeHandle& nh);
+
+ void publishMap();
+ void publishMapInflate(bool all_info = false);
+
+ void publishDepth();
+
+ bool hasDepthObservation();
+ bool odomValid();
+ void getRegion(Eigen::Vector3d& ori, Eigen::Vector3d& size);
+ inline double getResolution();
+ Eigen::Vector3d getOrigin();
+ int getVoxelNum();
+ bool getOdomDepthTimeout() { return md_.flag_depth_odom_timeout_; }
+
+ typedef std::shared_ptr Ptr;
+
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
+private:
+ MappingParameters mp_;
+ MappingData md_;
+
+ // get depth image and camera pose
+ void depthPoseCallback(const sensor_msgs::ImageConstPtr& img,
+ const geometry_msgs::PoseStampedConstPtr& pose);
+ void extrinsicCallback(const nav_msgs::OdometryConstPtr& odom);
+ void depthOdomCallback(const sensor_msgs::ImageConstPtr& img, const nav_msgs::OdometryConstPtr& odom);
+ void cloudCallback(const sensor_msgs::PointCloud2ConstPtr& img);
+ void odomCallback(const nav_msgs::OdometryConstPtr& odom);
+
+ // update occupancy by raycasting
+ void updateOccupancyCallback(const ros::TimerEvent& /*event*/);
+ void visCallback(const ros::TimerEvent& /*event*/);
+
+ // main update process
+ void projectDepthImage();
+ void raycastProcess();
+ void clearAndInflateLocalMap();
+
+ inline void inflatePoint(const Eigen::Vector3i& pt, int step, vector& pts);
+ int setCacheOccupancy(Eigen::Vector3d pos, int occ);
+ Eigen::Vector3d closetPointInMap(const Eigen::Vector3d& pt, const Eigen::Vector3d& camera_pt);
+
+ // typedef message_filters::sync_policies::ExactTime SyncPolicyImageOdom; typedef
+ // message_filters::sync_policies::ExactTime SyncPolicyImagePose;
+ typedef message_filters::sync_policies::ApproximateTime
+ SyncPolicyImageOdom;
+ typedef message_filters::sync_policies::ApproximateTime
+ SyncPolicyImagePose;
+ typedef shared_ptr> SynchronizerImagePose;
+ typedef shared_ptr> SynchronizerImageOdom;
+
+ ros::NodeHandle node_;
+ shared_ptr> depth_sub_;
+ shared_ptr> pose_sub_;
+ shared_ptr> odom_sub_;
+ SynchronizerImagePose sync_image_pose_;
+ SynchronizerImageOdom sync_image_odom_;
+
+ ros::Subscriber indep_cloud_sub_, indep_odom_sub_, extrinsic_sub_;
+ ros::Publisher map_pub_, map_inf_pub_;
+ ros::Timer occ_timer_, vis_timer_;
+
+ //
+ uniform_real_distribution rand_noise_;
+ normal_distribution rand_noise2_;
+ default_random_engine eng_;
+};
+
+/* ============================== definition of inline function
+ * ============================== */
+
+inline int GridMap::toAddress(const Eigen::Vector3i& id) {
+ return id(0) * mp_.map_voxel_num_(1) * mp_.map_voxel_num_(2) + id(1) * mp_.map_voxel_num_(2) + id(2);
+}
+
+inline int GridMap::toAddress(int& x, int& y, int& z) {
+ return x * mp_.map_voxel_num_(1) * mp_.map_voxel_num_(2) + y * mp_.map_voxel_num_(2) + z;
+}
+
+inline void GridMap::boundIndex(Eigen::Vector3i& id) {
+ Eigen::Vector3i id1;
+ id1(0) = max(min(id(0), mp_.map_voxel_num_(0) - 1), 0);
+ id1(1) = max(min(id(1), mp_.map_voxel_num_(1) - 1), 0);
+ id1(2) = max(min(id(2), mp_.map_voxel_num_(2) - 1), 0);
+ id = id1;
+}
+
+inline bool GridMap::isUnknown(const Eigen::Vector3i& id) {
+ Eigen::Vector3i id1 = id;
+ boundIndex(id1);
+ return md_.occupancy_buffer_[toAddress(id1)] < mp_.clamp_min_log_ - 1e-3;
+}
+
+inline bool GridMap::isUnknown(const Eigen::Vector3d& pos) {
+ Eigen::Vector3i idc;
+ posToIndex(pos, idc);
+ return isUnknown(idc);
+}
+
+inline bool GridMap::isKnownFree(const Eigen::Vector3i& id) {
+ Eigen::Vector3i id1 = id;
+ boundIndex(id1);
+ int adr = toAddress(id1);
+
+ // return md_.occupancy_buffer_[adr] >= mp_.clamp_min_log_ &&
+ // md_.occupancy_buffer_[adr] < mp_.min_occupancy_log_;
+ return md_.occupancy_buffer_[adr] >= mp_.clamp_min_log_ && md_.occupancy_buffer_inflate_[adr] == 0;
+}
+
+inline bool GridMap::isKnownOccupied(const Eigen::Vector3i& id) {
+ Eigen::Vector3i id1 = id;
+ boundIndex(id1);
+ int adr = toAddress(id1);
+
+ return md_.occupancy_buffer_inflate_[adr] == 1;
+}
+
+inline void GridMap::setOccupied(Eigen::Vector3d pos) {
+ if (!isInMap(pos)) return;
+
+ Eigen::Vector3i id;
+ posToIndex(pos, id);
+
+ md_.occupancy_buffer_inflate_[id(0) * mp_.map_voxel_num_(1) * mp_.map_voxel_num_(2) +
+ id(1) * mp_.map_voxel_num_(2) + id(2)] = 1;
+}
+
+inline void GridMap::setOccupancy(Eigen::Vector3d pos, double occ) {
+ if (occ != 1 && occ != 0) {
+ cout << "occ value error!" << endl;
+ return;
+ }
+
+ if (!isInMap(pos)) return;
+
+ Eigen::Vector3i id;
+ posToIndex(pos, id);
+
+ md_.occupancy_buffer_[toAddress(id)] = occ;
+}
+
+inline int GridMap::getOccupancy(Eigen::Vector3d pos) {
+ if (!isInMap(pos)) return -1;
+
+ Eigen::Vector3i id;
+ posToIndex(pos, id);
+
+ return md_.occupancy_buffer_[toAddress(id)] > mp_.min_occupancy_log_ ? 1 : 0;
+}
+
+inline int GridMap::getInflateOccupancy(Eigen::Vector3d pos) {
+ if (!isInMap(pos)) return -1;
+
+ Eigen::Vector3i id;
+ posToIndex(pos, id);
+
+ return int(md_.occupancy_buffer_inflate_[toAddress(id)]);
+}
+
+inline int GridMap::getOccupancy(Eigen::Vector3i id) {
+ if (id(0) < 0 || id(0) >= mp_.map_voxel_num_(0) || id(1) < 0 || id(1) >= mp_.map_voxel_num_(1) ||
+ id(2) < 0 || id(2) >= mp_.map_voxel_num_(2))
+ return -1;
+
+ return md_.occupancy_buffer_[toAddress(id)] > mp_.min_occupancy_log_ ? 1 : 0;
+}
+
+inline bool GridMap::isInMap(const Eigen::Vector3d& pos) {
+ if (pos(0) < mp_.map_min_boundary_(0) + 1e-4 || pos(1) < mp_.map_min_boundary_(1) + 1e-4 ||
+ pos(2) < mp_.map_min_boundary_(2) + 1e-4) {
+ // cout << "less than min range!" << endl;
+ return false;
+ }
+ if (pos(0) > mp_.map_max_boundary_(0) - 1e-4 || pos(1) > mp_.map_max_boundary_(1) - 1e-4 ||
+ pos(2) > mp_.map_max_boundary_(2) - 1e-4) {
+ return false;
+ }
+ return true;
+}
+
+inline bool GridMap::isInMap(const Eigen::Vector3i& idx) {
+ if (idx(0) < 0 || idx(1) < 0 || idx(2) < 0) {
+ return false;
+ }
+ if (idx(0) > mp_.map_voxel_num_(0) - 1 || idx(1) > mp_.map_voxel_num_(1) - 1 ||
+ idx(2) > mp_.map_voxel_num_(2) - 1) {
+ return false;
+ }
+ return true;
+}
+
+inline void GridMap::posToIndex(const Eigen::Vector3d& pos, Eigen::Vector3i& id) {
+ for (int i = 0; i < 3; ++i) id(i) = floor((pos(i) - mp_.map_origin_(i)) * mp_.resolution_inv_);
+}
+
+inline void GridMap::indexToPos(const Eigen::Vector3i& id, Eigen::Vector3d& pos) {
+ for (int i = 0; i < 3; ++i) pos(i) = (id(i) + 0.5) * mp_.resolution_ + mp_.map_origin_(i);
+}
+
+inline void GridMap::inflatePoint(const Eigen::Vector3i& pt, int step, vector& pts) {
+ int num = 0;
+
+ /* ---------- + shape inflate ---------- */
+ // for (int x = -step; x <= step; ++x)
+ // {
+ // if (x == 0)
+ // continue;
+ // pts[num++] = Eigen::Vector3i(pt(0) + x, pt(1), pt(2));
+ // }
+ // for (int y = -step; y <= step; ++y)
+ // {
+ // if (y == 0)
+ // continue;
+ // pts[num++] = Eigen::Vector3i(pt(0), pt(1) + y, pt(2));
+ // }
+ // for (int z = -1; z <= 1; ++z)
+ // {
+ // pts[num++] = Eigen::Vector3i(pt(0), pt(1), pt(2) + z);
+ // }
+
+ /* ---------- all inflate ---------- */
+ for (int x = -step; x <= step; ++x)
+ for (int y = -step; y <= step; ++y)
+ for (int z = -step; z <= step; ++z) {
+ pts[num++] = Eigen::Vector3i(pt(0) + x, pt(1) + y, pt(2) + z);
+ }
+}
+
+inline double GridMap::getResolution() { return mp_.resolution_; }
+
+#endif
+
+// #ifndef _GRID_MAP_H
+// #define _GRID_MAP_H
+
+// #include
+// #include
+// #include
+// #include
+// #include
+// #include
+// #include
+// #include
+// #include
+// #include
+// #include
+
+// #include
+// #include
+// #include
+
+// #include
+// #include
+// #include
+// #include
+
+// #include
+
+// #define logit(x) (log((x) / (1 - (x))))
+
+// using namespace std;
+
+// // voxel hashing
+// template
+// struct matrix_hash : std::unary_function {
+// std::size_t operator()(T const& matrix) const {
+// size_t seed = 0;
+// for (size_t i = 0; i < matrix.size(); ++i) {
+// auto elem = *(matrix.data() + i);
+// seed ^= std::hash()(elem) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+// }
+// return seed;
+// }
+// };
+
+// // constant parameters
+
+// struct MappingParameters {
+
+// /* map properties */
+// Eigen::Vector3d map_origin_, map_size_;
+// Eigen::Vector3d map_min_boundary_, map_max_boundary_; // map range in pos
+// Eigen::Vector3i map_voxel_num_; // map range in index
+// Eigen::Vector3d local_update_range_;
+// double resolution_, resolution_inv_;
+// double obstacles_inflation_;
+// string frame_id_;
+// int pose_type_;
+
+// /* camera parameters */
+// double cx_, cy_, fx_, fy_;
+
+// /* depth image projection filtering */
+// double depth_filter_maxdist_, depth_filter_mindist_, depth_filter_tolerance_;
+// int depth_filter_margin_;
+// bool use_depth_filter_;
+// double k_depth_scaling_factor_;
+// int skip_pixel_;
+
+// /* raycasting */
+// double p_hit_, p_miss_, p_min_, p_max_, p_occ_; // occupancy probability
+// double prob_hit_log_, prob_miss_log_, clamp_min_log_, clamp_max_log_,
+// min_occupancy_log_; // logit of occupancy probability
+// double min_ray_length_, max_ray_length_; // range of doing raycasting
+
+// /* local map update and clear */
+// int local_map_margin_;
+
+// /* visualization and computation time display */
+// double visualization_truncate_height_, virtual_ceil_height_, ground_height_;
+// bool show_occ_time_;
+
+// /* active mapping */
+// double unknown_flag_;
+// };
+
+// // intermediate mapping data for fusion
+
+// struct MappingData {
+// // main map data, occupancy of each voxel and Euclidean distance
+
+// std::vector occupancy_buffer_;
+// std::vector occupancy_buffer_inflate_;
+
+// // camera position and pose data
+
+// Eigen::Vector3d camera_pos_, last_camera_pos_;
+// Eigen::Quaterniond camera_q_, last_camera_q_;
+
+// // depth image data
+
+// cv::Mat depth_image_, last_depth_image_;
+// int image_cnt_;
+
+// // flags of map state
+
+// bool occ_need_update_, local_updated_;
+// bool has_first_depth_;
+// bool has_odom_, has_cloud_;
+
+// // depth image projected point cloud
+
+// vector proj_points_;
+// int proj_points_cnt;
+
+// // flag buffers for speeding up raycasting
+
+// vector count_hit_, count_hit_and_miss_;
+// vector flag_traverse_, flag_rayend_;
+// char raycast_num_;
+// queue cache_voxel_;
+
+// // range of updating grid
+
+// Eigen::Vector3i local_bound_min_, local_bound_max_;
+
+// // computation time
+
+// double fuse_time_, max_fuse_time_;
+// int update_num_;
+
+// EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+// };
+
+// class GridMap {
+// public:
+// GridMap() {}
+// ~GridMap() {}
+
+// enum { POSE_STAMPED = 1, ODOMETRY = 2, INVALID_IDX = -10000 };
+
+// // occupancy map management
+// void resetBuffer();
+// void resetBuffer(Eigen::Vector3d min, Eigen::Vector3d max);
+
+// inline void posToIndex(const Eigen::Vector3d& pos, Eigen::Vector3i& id);
+// inline void indexToPos(const Eigen::Vector3i& id, Eigen::Vector3d& pos);
+// inline int toAddress(const Eigen::Vector3i& id);
+// inline int toAddress(int& x, int& y, int& z);
+// inline bool isInMap(const Eigen::Vector3d& pos);
+// inline bool isInMap(const Eigen::Vector3i& idx);
+
+// inline void setOccupancy(Eigen::Vector3d pos, double occ = 1);
+// inline void setOccupied(Eigen::Vector3d pos);
+// inline int getOccupancy(Eigen::Vector3d pos);
+// inline int getOccupancy(Eigen::Vector3i id);
+// inline int getInflateOccupancy(Eigen::Vector3d pos);
+
+// inline void boundIndex(Eigen::Vector3i& id);
+// inline bool isUnknown(const Eigen::Vector3i& id);
+// inline bool isUnknown(const Eigen::Vector3d& pos);
+// inline bool isKnownFree(const Eigen::Vector3i& id);
+// inline bool isKnownOccupied(const Eigen::Vector3i& id);
+
+// void initMap(ros::NodeHandle& nh);
+
+// void publishMap();
+// void publishMapInflate(bool all_info = false);
+
+// void publishUnknown();
+// void publishDepth();
+
+// bool hasDepthObservation();
+// bool odomValid();
+// void getRegion(Eigen::Vector3d& ori, Eigen::Vector3d& size);
+// inline double getResolution();
+// Eigen::Vector3d getOrigin();
+// int getVoxelNum();
+
+// typedef std::shared_ptr Ptr;
+
+// EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
+// private:
+// MappingParameters mp_;
+// MappingData md_;
+
+// // get depth image and camera pose
+// void depthPoseCallback(const sensor_msgs::ImageConstPtr& img,
+// const geometry_msgs::PoseStampedConstPtr& pose);
+// void depthOdomCallback(const sensor_msgs::ImageConstPtr& img, const nav_msgs::OdometryConstPtr& odom);
+// void cloudCallback(const sensor_msgs::PointCloud2ConstPtr& img);
+// void odomCallback(const nav_msgs::OdometryConstPtr& odom);
+
+// // update occupancy by raycasting
+// void updateOccupancyCallback(const ros::TimerEvent& /*event*/);
+// void visCallback(const ros::TimerEvent& /*event*/);
+
+// // main update process
+// void projectDepthImage();
+// void raycastProcess();
+// void clearAndInflateLocalMap();
+
+// inline void inflatePoint(const Eigen::Vector3i& pt, int step, vector& pts);
+// int setCacheOccupancy(Eigen::Vector3d pos, int occ);
+// Eigen::Vector3d closetPointInMap(const Eigen::Vector3d& pt, const Eigen::Vector3d& camera_pt);
+
+// // typedef message_filters::sync_policies::ExactTime SyncPolicyImageOdom; typedef
+// // message_filters::sync_policies::ExactTime SyncPolicyImagePose;
+// typedef message_filters::sync_policies::ApproximateTime
+// SyncPolicyImageOdom;
+// typedef message_filters::sync_policies::ApproximateTime
+// SyncPolicyImagePose;
+// typedef shared_ptr> SynchronizerImagePose;
+// typedef shared_ptr> SynchronizerImageOdom;
+
+// ros::NodeHandle node_;
+// shared_ptr> depth_sub_;
+// shared_ptr> pose_sub_;
+// shared_ptr> odom_sub_;
+// SynchronizerImagePose sync_image_pose_;
+// SynchronizerImageOdom sync_image_odom_;
+
+// ros::Subscriber indep_cloud_sub_, indep_odom_sub_;
+// ros::Publisher map_pub_, map_inf_pub_;
+// ros::Publisher unknown_pub_;
+// ros::Timer occ_timer_, vis_timer_;
+
+// //
+// uniform_real_distribution rand_noise_;
+// normal_distribution rand_noise2_;
+// default_random_engine eng_;
+// };
+
+// /* ============================== definition of inline function
+// * ============================== */
+
+// inline int GridMap::toAddress(const Eigen::Vector3i& id) {
+// return id(0) * mp_.map_voxel_num_(1) * mp_.map_voxel_num_(2) + id(1) * mp_.map_voxel_num_(2) + id(2);
+// }
+
+// inline int GridMap::toAddress(int& x, int& y, int& z) {
+// return x * mp_.map_voxel_num_(1) * mp_.map_voxel_num_(2) + y * mp_.map_voxel_num_(2) + z;
+// }
+
+// inline void GridMap::boundIndex(Eigen::Vector3i& id) {
+// Eigen::Vector3i id1;
+// id1(0) = max(min(id(0), mp_.map_voxel_num_(0) - 1), 0);
+// id1(1) = max(min(id(1), mp_.map_voxel_num_(1) - 1), 0);
+// id1(2) = max(min(id(2), mp_.map_voxel_num_(2) - 1), 0);
+// id = id1;
+// }
+
+// inline bool GridMap::isUnknown(const Eigen::Vector3i& id) {
+// Eigen::Vector3i id1 = id;
+// boundIndex(id1);
+// return md_.occupancy_buffer_[toAddress(id1)] < mp_.clamp_min_log_ - 1e-3;
+// }
+
+// inline bool GridMap::isUnknown(const Eigen::Vector3d& pos) {
+// Eigen::Vector3i idc;
+// posToIndex(pos, idc);
+// return isUnknown(idc);
+// }
+
+// inline bool GridMap::isKnownFree(const Eigen::Vector3i& id) {
+// Eigen::Vector3i id1 = id;
+// boundIndex(id1);
+// int adr = toAddress(id1);
+
+// // return md_.occupancy_buffer_[adr] >= mp_.clamp_min_log_ &&
+// // md_.occupancy_buffer_[adr] < mp_.min_occupancy_log_;
+// return md_.occupancy_buffer_[adr] >= mp_.clamp_min_log_ && md_.occupancy_buffer_inflate_[adr] == 0;
+// }
+
+// inline bool GridMap::isKnownOccupied(const Eigen::Vector3i& id) {
+// Eigen::Vector3i id1 = id;
+// boundIndex(id1);
+// int adr = toAddress(id1);
+
+// return md_.occupancy_buffer_inflate_[adr] == 1;
+// }
+
+// inline void GridMap::setOccupied(Eigen::Vector3d pos) {
+// if (!isInMap(pos)) return;
+
+// Eigen::Vector3i id;
+// posToIndex(pos, id);
+
+// md_.occupancy_buffer_inflate_[id(0) * mp_.map_voxel_num_(1) * mp_.map_voxel_num_(2) +
+// id(1) * mp_.map_voxel_num_(2) + id(2)] = 1;
+// }
+
+// inline void GridMap::setOccupancy(Eigen::Vector3d pos, double occ) {
+// if (occ != 1 && occ != 0) {
+// cout << "occ value error!" << endl;
+// return;
+// }
+
+// if (!isInMap(pos)) return;
+
+// Eigen::Vector3i id;
+// posToIndex(pos, id);
+
+// md_.occupancy_buffer_[toAddress(id)] = occ;
+// }
+
+// inline int GridMap::getOccupancy(Eigen::Vector3d pos) {
+// if (!isInMap(pos)) return -1;
+
+// Eigen::Vector3i id;
+// posToIndex(pos, id);
+
+// return md_.occupancy_buffer_[toAddress(id)] > mp_.min_occupancy_log_ ? 1 : 0;
+// }
+
+// inline int GridMap::getInflateOccupancy(Eigen::Vector3d pos) {
+// if (!isInMap(pos)) return -1;
+
+// Eigen::Vector3i id;
+// posToIndex(pos, id);
+
+// return int(md_.occupancy_buffer_inflate_[toAddress(id)]);
+// }
+
+// inline int GridMap::getOccupancy(Eigen::Vector3i id) {
+// if (id(0) < 0 || id(0) >= mp_.map_voxel_num_(0) || id(1) < 0 || id(1) >= mp_.map_voxel_num_(1) ||
+// id(2) < 0 || id(2) >= mp_.map_voxel_num_(2))
+// return -1;
+
+// return md_.occupancy_buffer_[toAddress(id)] > mp_.min_occupancy_log_ ? 1 : 0;
+// }
+
+// inline bool GridMap::isInMap(const Eigen::Vector3d& pos) {
+// if (pos(0) < mp_.map_min_boundary_(0) + 1e-4 || pos(1) < mp_.map_min_boundary_(1) + 1e-4 ||
+// pos(2) < mp_.map_min_boundary_(2) + 1e-4) {
+// // cout << "less than min range!" << endl;
+// return false;
+// }
+// if (pos(0) > mp_.map_max_boundary_(0) - 1e-4 || pos(1) > mp_.map_max_boundary_(1) - 1e-4 ||
+// pos(2) > mp_.map_max_boundary_(2) - 1e-4) {
+// return false;
+// }
+// return true;
+// }
+
+// inline bool GridMap::isInMap(const Eigen::Vector3i& idx) {
+// if (idx(0) < 0 || idx(1) < 0 || idx(2) < 0) {
+// return false;
+// }
+// if (idx(0) > mp_.map_voxel_num_(0) - 1 || idx(1) > mp_.map_voxel_num_(1) - 1 ||
+// idx(2) > mp_.map_voxel_num_(2) - 1) {
+// return false;
+// }
+// return true;
+// }
+
+// inline void GridMap::posToIndex(const Eigen::Vector3d& pos, Eigen::Vector3i& id) {
+// for (int i = 0; i < 3; ++i) id(i) = floor((pos(i) - mp_.map_origin_(i)) * mp_.resolution_inv_);
+// }
+
+// inline void GridMap::indexToPos(const Eigen::Vector3i& id, Eigen::Vector3d& pos) {
+// for (int i = 0; i < 3; ++i) pos(i) = (id(i) + 0.5) * mp_.resolution_ + mp_.map_origin_(i);
+// }
+
+// inline void GridMap::inflatePoint(const Eigen::Vector3i& pt, int step, vector& pts) {
+// int num = 0;
+
+// /* ---------- + shape inflate ---------- */
+// // for (int x = -step; x <= step; ++x)
+// // {
+// // if (x == 0)
+// // continue;
+// // pts[num++] = Eigen::Vector3i(pt(0) + x, pt(1), pt(2));
+// // }
+// // for (int y = -step; y <= step; ++y)
+// // {
+// // if (y == 0)
+// // continue;
+// // pts[num++] = Eigen::Vector3i(pt(0), pt(1) + y, pt(2));
+// // }
+// // for (int z = -1; z <= 1; ++z)
+// // {
+// // pts[num++] = Eigen::Vector3i(pt(0), pt(1), pt(2) + z);
+// // }
+
+// /* ---------- all inflate ---------- */
+// for (int x = -step; x <= step; ++x)
+// for (int y = -step; y <= step; ++y)
+// for (int z = -step; z <= step; ++z) {
+// pts[num++] = Eigen::Vector3i(pt(0) + x, pt(1) + y, pt(2) + z);
+// }
+// }
+
+// inline double GridMap::getResolution() { return mp_.resolution_; }
+
+// #endif
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/linear_obj_model.hpp b/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/linear_obj_model.hpp
new file mode 100644
index 0000000..574e3b0
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/linear_obj_model.hpp
@@ -0,0 +1,267 @@
+/**
+* This file is part of Fast-Planner.
+*
+* Copyright 2019 Boyu Zhou, Aerial Robotics Group, Hong Kong University of Science and Technology,
+* Developed by Boyu Zhou ,
+* for more information see .
+* If you use this code, please cite the respective publications as
+* listed on the above website.
+*
+* Fast-Planner is free software: you can redistribute it and/or modify
+* it under the terms of the GNU Lesser General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* Fast-Planner is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with Fast-Planner. If not, see .
+*/
+
+
+
+#ifndef _LINEAR_OBJ_MODEL_H_
+#define _LINEAR_OBJ_MODEL_H_
+
+#include
+#include
+#include
+
+class LinearObjModel {
+private:
+ bool last_out_bound_{false};
+ int input_type_;
+ /* data */
+public:
+ LinearObjModel(/* args */);
+ ~LinearObjModel();
+
+ void initialize(Eigen::Vector3d p, Eigen::Vector3d v, Eigen::Vector3d a, double yaw, double yaw_dot,
+ Eigen::Vector3d color, Eigen::Vector3d scale, int input_type);
+
+ void setLimits(Eigen::Vector3d bound, Eigen::Vector2d vel, Eigen::Vector2d acc);
+
+ void update(double dt); // linear trippler integrator model
+
+ static bool collide(LinearObjModel& obj1, LinearObjModel& obj2);
+
+ // void setInput(Eigen::Vector3d acc) {
+ // acc_ = acc;
+ // }
+ void setInput(Eigen::Vector3d vel) {
+ vel_ = vel;
+ }
+
+ void setYawDot(double yaw_dot) {
+ yaw_dot_ = yaw_dot;
+ }
+
+ Eigen::Vector3d getPosition() {
+ return pos_;
+ }
+ void setPosition(Eigen::Vector3d pos) {
+ pos_ = pos;
+ }
+
+ Eigen::Vector3d getVelocity() {
+ return vel_;
+ }
+
+ void setVelocity(double x, double y, double z) {
+ vel_ = Eigen::Vector3d(x, y, z);
+ }
+
+ Eigen::Vector3d getColor() {
+ return color_;
+ }
+ Eigen::Vector3d getScale() {
+ return scale_;
+ }
+
+ double getYaw() {
+ return yaw_;
+ }
+
+private:
+ Eigen::Vector3d pos_, vel_, acc_;
+ Eigen::Vector3d color_, scale_;
+ double yaw_, yaw_dot_;
+
+ Eigen::Vector3d bound_;
+ Eigen::Vector2d limit_v_, limit_a_;
+};
+
+LinearObjModel::LinearObjModel(/* args */) {
+}
+
+LinearObjModel::~LinearObjModel() {
+}
+
+void LinearObjModel::initialize(Eigen::Vector3d p, Eigen::Vector3d v, Eigen::Vector3d a, double yaw,
+ double yaw_dot, Eigen::Vector3d color, Eigen::Vector3d scale, int input_type) {
+ pos_ = p;
+ vel_ = v;
+ acc_ = a;
+ color_ = color;
+ scale_ = scale;
+ input_type_ = input_type;
+
+ yaw_ = yaw;
+ yaw_dot_ = yaw_dot;
+}
+
+void LinearObjModel::setLimits(Eigen::Vector3d bound, Eigen::Vector2d vel, Eigen::Vector2d acc) {
+ bound_ = bound;
+ limit_v_ = vel;
+ limit_a_ = acc;
+}
+
+void LinearObjModel::update(double dt) {
+ Eigen::Vector3d p0, v0, a0;
+ p0 = pos_, v0 = vel_, a0 = acc_;
+ //std::cout << v0.transpose() << std::endl;
+
+ /* ---------- use acc as input ---------- */
+ if ( input_type_ == 2 )
+ {
+ vel_ = v0 + acc_ * dt;
+ for (int i = 0; i < 3; ++i)
+ {
+ if (vel_(i) > 0) vel_(i) = std::max(limit_v_(0), std::min(vel_(i),
+ limit_v_(1)));
+ if (vel_(i) <= 0) vel_(i) = std::max(-limit_v_(1), std::min(vel_(i),
+ -limit_v_(0)));
+ }
+
+ pos_ = p0 + v0 * dt + 0.5 * acc_ * pow(dt, 2);
+
+ /* ---------- reflect acc when collide with bound ---------- */
+ if ( pos_(0) <= bound_(0) && pos_(0) >= -bound_(0) &&
+ pos_(1) <= bound_(1) && pos_(1) >= -bound_(1) &&
+ pos_(2) <= bound_(2) && pos_(2) >= 0
+ )
+ {
+ last_out_bound_ = false;
+ }
+ else if ( !last_out_bound_ )
+ {
+ last_out_bound_ = true;
+
+ // if ( pos_(0) > bound_(0) || pos_(0) < -bound_(0) ) acc_(0) = -acc_(0);
+ // if ( pos_(1) > bound_(1) || pos_(1) < -bound_(1) ) acc_(1) = -acc_(1);
+ // if ( pos_(2) > bound_(2) || pos_(2) < -bound_(2) ) acc_(2) = -acc_(2);
+ acc_ = -acc_;
+ //ROS_ERROR("AAAAAAAAAAAAAAAAAAa");
+ }
+ }
+ // for (int i = 0; i < 2; ++i)
+ // {
+ // pos_(i) = std::min(pos_(i), bound_(i));
+ // pos_(i) = std::max(pos_(i), -bound_(i));
+ // }
+ // pos_(2) = std::min(pos_(2), bound_(2));
+ // pos_(2) = std::max(pos_(2), 0.0);
+
+ /* ---------- use vel as input ---------- */
+ else if ( input_type_ == 1 )
+ {
+ pos_ = p0 + v0 * dt;
+ for (int i = 0; i < 2; ++i) {
+ pos_(i) = std::min(pos_(i), bound_(i));
+ pos_(i) = std::max(pos_(i), -bound_(i));
+ }
+ pos_(2) = std::min(pos_(2), bound_(2));
+ pos_(2) = std::max(pos_(2), 0.0);
+
+ yaw_ += yaw_dot_ * dt;
+
+ const double PI = 3.1415926;
+ if (yaw_ > 2 * PI) yaw_ -= 2 * PI;
+
+ const double tol = 0.1;
+ if (pos_(0) > bound_(0) - tol) {
+ pos_(0) = bound_(0) - tol;
+ vel_(0) = -vel_(0);
+ }
+ if (pos_(0) < -bound_(0) + tol) {
+ pos_(0) = -bound_(0) + tol;
+ vel_(0) = -vel_(0);
+ }
+
+ if (pos_(1) > bound_(1) - tol) {
+ pos_(1) = bound_(1) - tol;
+ vel_(1) = -vel_(1);
+ }
+ if (pos_(1) < -bound_(1) + tol) {
+ pos_(1) = -bound_(1) + tol;
+ vel_(1) = -vel_(1);
+ }
+
+ if (pos_(2) > bound_(2) - tol) {
+ pos_(2) = bound_(2) - tol;
+ vel_(2) = -vel_(2);
+ }
+ if (pos_(2) < tol) {
+ pos_(2) = tol;
+ vel_(2) = -vel_(2);
+ }
+ }
+
+ // /* ---------- reflect when collide with bound ---------- */
+
+
+ //std::cout << pos_.transpose() << " " << bound_.transpose() << std::endl;
+}
+
+bool LinearObjModel::collide(LinearObjModel& obj1, LinearObjModel& obj2) {
+ Eigen::Vector3d pos1, pos2, vel1, vel2, scale1, scale2;
+ pos1 = obj1.getPosition();
+ vel1 = obj1.getVelocity();
+ scale1 = obj1.getScale();
+
+ pos2 = obj2.getPosition();
+ vel2 = obj2.getVelocity();
+ scale2 = obj2.getScale();
+
+ /* ---------- collide ---------- */
+ bool collide = fabs(pos1(0) - pos2(0)) < 0.5 * (scale1(0) + scale2(0)) &&
+ fabs(pos1(1) - pos2(1)) < 0.5 * (scale1(1) + scale2(1)) &&
+ fabs(pos1(2) - pos2(2)) < 0.5 * (scale1(2) + scale2(2));
+
+ if (collide) {
+ double tol[3];
+ tol[0] = 0.5 * (scale1(0) + scale2(0)) - fabs(pos1(0) - pos2(0));
+ tol[1] = 0.5 * (scale1(1) + scale2(1)) - fabs(pos1(1) - pos2(1));
+ tol[2] = 0.5 * (scale1(2) + scale2(2)) - fabs(pos1(2) - pos2(2));
+
+ for (int i = 0; i < 3; ++i) {
+ if (tol[i] < tol[(i + 1) % 3] && tol[i] < tol[(i + 2) % 3]) {
+ vel1(i) = -vel1(i);
+ vel2(i) = -vel2(i);
+ obj1.setVelocity(vel1(0), vel1(1), vel1(2));
+ obj2.setVelocity(vel2(0), vel2(1), vel2(2));
+
+ if (pos1(i) >= pos2(i)) {
+ pos1(i) += tol[i];
+ pos2(i) -= tol[i];
+ } else {
+ pos1(i) -= tol[i];
+ pos2(i) += tol[i];
+ }
+ obj1.setPosition(pos1);
+ obj2.setPosition(pos2);
+
+ break;
+ }
+ }
+
+ return true;
+ } else {
+ return false;
+ }
+}
+
+#endif
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/obj_predictor.h b/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/obj_predictor.h
new file mode 100644
index 0000000..fb2e86b
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/obj_predictor.h
@@ -0,0 +1,176 @@
+/**
+* This file is part of Fast-Planner.
+*
+* Copyright 2019 Boyu Zhou, Aerial Robotics Group, Hong Kong University of Science and Technology,
+* Developed by Boyu Zhou ,
+* for more information see .
+* If you use this code, please cite the respective publications as
+* listed on the above website.
+*
+* Fast-Planner is free software: you can redistribute it and/or modify
+* it under the terms of the GNU Lesser General Public License as published by
+* the Free Software Foundation, either version 3 of the License, or
+* (at your option) any later version.
+*
+* Fast-Planner is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+* GNU General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with Fast-Planner. If not, see .
+*/
+
+
+
+#ifndef _OBJ_PREDICTOR_H_
+#define _OBJ_PREDICTOR_H_
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+using std::cout;
+using std::endl;
+using std::list;
+using std::shared_ptr;
+using std::unique_ptr;
+using std::vector;
+
+namespace fast_planner {
+class PolynomialPrediction;
+typedef shared_ptr> ObjPrediction;
+typedef shared_ptr> ObjScale;
+
+/* ========== prediction polynomial ========== */
+class PolynomialPrediction {
+private:
+ vector> polys;
+ double t1, t2; // start / end
+ ros::Time global_start_time_;
+
+public:
+ PolynomialPrediction(/* args */) {
+ }
+ ~PolynomialPrediction() {
+ }
+
+ void setPolynomial(vector>& pls) {
+ polys = pls;
+ }
+ void setTime(double t1, double t2) {
+ this->t1 = t1;
+ this->t2 = t2;
+ }
+ void setGlobalStartTime(ros::Time global_start_time) {
+ global_start_time_ = global_start_time;
+ }
+
+ bool valid() {
+ return polys.size() == 3;
+ }
+
+ /* note that t should be in [t1, t2] */
+ Eigen::Vector3d evaluate(double t) {
+ Eigen::Matrix tv;
+ tv << 1.0, pow(t, 1), pow(t, 2), pow(t, 3), pow(t, 4), pow(t, 5);
+
+ Eigen::Vector3d pt;
+ pt(0) = tv.dot(polys[0]), pt(1) = tv.dot(polys[1]), pt(2) = tv.dot(polys[2]);
+
+ return pt;
+ }
+
+ Eigen::Vector3d evaluateConstVel(double t) {
+ Eigen::Matrix tv;
+ tv << 1.0, pow(t-global_start_time_.toSec(), 1);
+
+ // cout << t-global_start_time_.toSec() << endl;
+
+ Eigen::Vector3d pt;
+ pt(0) = tv.dot(polys[0].head(2)), pt(1) = tv.dot(polys[1].head(2)), pt(2) = tv.dot(polys[2].head(2));
+
+ return pt;
+ }
+};
+
+/* ========== subscribe and record object history ========== */
+class ObjHistory {
+public:
+ int skip_num_;
+ int queue_size_;
+ ros::Time global_start_time_;
+
+ ObjHistory() {
+ }
+ ~ObjHistory() {
+ }
+
+ void init(int id, int skip_num, int queue_size, ros::Time global_start_time);
+
+ void poseCallback(const geometry_msgs::PoseStampedConstPtr& msg);
+
+ void clear() {
+ history_.clear();
+ }
+
+ void getHistory(list& his) {
+ his = history_;
+ }
+
+private:
+ list history_; // x,y,z;t
+ int skip_;
+ int obj_idx_;
+ Eigen::Vector3d scale_;
+};
+
+/* ========== predict future trajectory using history ========== */
+class ObjPredictor {
+private:
+ ros::NodeHandle node_handle_;
+
+ int obj_num_;
+ double lambda_;
+ double predict_rate_;
+
+ vector pose_subs_;
+ ros::Subscriber marker_sub_;
+ ros::Timer predict_timer_;
+ vector> obj_histories_;
+
+ /* share data with planner */
+ ObjPrediction predict_trajs_;
+ ObjScale obj_scale_;
+ vector scale_init_;
+
+ void markerCallback(const visualization_msgs::MarkerConstPtr& msg);
+
+ void predictCallback(const ros::TimerEvent& e);
+ void predictPolyFit();
+ void predictConstVel();
+
+public:
+ ObjPredictor(/* args */);
+ ObjPredictor(ros::NodeHandle& node);
+ ~ObjPredictor();
+
+ void init();
+
+ ObjPrediction getPredictionTraj();
+ ObjScale getObjScale();
+ int getObjNums() {return obj_num_;}
+
+ Eigen::Vector3d evaluatePoly(int obs_id, double time);
+ Eigen::Vector3d evaluateConstVel(int obs_id, double time);
+
+ typedef shared_ptr Ptr;
+};
+
+} // namespace fast_planner
+
+#endif
\ No newline at end of file
diff --git a/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/raycast.h b/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/raycast.h
new file mode 100644
index 0000000..23e0876
--- /dev/null
+++ b/motion_planning/3d/ego_planner/planner/plan_env/include/plan_env/raycast.h
@@ -0,0 +1,63 @@
+#ifndef RAYCAST_H_
+#define RAYCAST_H_
+
+#include