25#include <Eigen/Geometry>
46 core::GlobalTrajectory gtraj;
47 for (
const auto& p : traj)
49 core::GlobalTrajectoryPoint gp;
51 Eigen::Translation3f{Eigen::Vector3f(
static_cast<float>(p.x),
52 static_cast<float>(p.y),
54 Eigen::AngleAxisf{
static_cast<float>(p.theta), Eigen::Vector3f::UnitZ()};
55 gp.velocity =
static_cast<float>(p.v);
56 gtraj.mutablePoints().push_back(gp);
63 std::vector<CenterPoint> traj;
64 for (
const auto& gp : gtraj.points())
67 p.
x = gp.waypoint.pose.translation().x();
68 p.y = gp.waypoint.pose.translation().y();
69 const auto& R = gp.waypoint.pose.rotation();
70 p.theta = std::atan2(R(1, 0), R(0, 0));
80 wrapper(&x, &y, &out);
84 bool isInsideBounds(
double x,
double y,
const algorithms::Costmap& costmap)
86 const auto local = costmap.origin().inverse() * Eigen::Vector2f{
static_cast<float>(x),
87 static_cast<float>(y)};
88 const auto& bounds = costmap.getLocalSceneBounds();
89 return local.x() >= bounds.min.x() && local.x() <= bounds.max.x() &&
90 local.y() >= bounds.min.y() && local.y() <= bounds.max.y();
93 void pushWaypointsToSafeZone(std::vector<CenterPoint>& traj,
95 const algorithms::Costmap& costmap,
96 const std::vector<CenterPoint>& originalTargets,
102 const std::vector<std::pair<double, double>> neighbours = {
108 {stepSize, stepSize},
109 {stepSize, -stepSize},
110 {-stepSize, stepSize},
111 {-stepSize, -stepSize},
114 for (
int idx = 1; idx + 1 <
static_cast<int>(traj.size()); ++idx)
117 const double origX = originalTargets[idx].x;
118 const double origY = originalTargets[idx].y;
120 auto withinDeviation = [&](
double x,
double y)
122 return std::hypot(x - origX, y - origY) <= maxDeviation + 1e-6;
125 double bestDist = queryDistance(wrapper, p.x, p.y);
127 if (bestDist >= clearance)
130 for (
int iter = 0; iter < maxIterations; ++iter)
134 bool improved =
false;
136 for (
const auto& [dx, dy] : neighbours)
138 double candX = p.x + dx;
139 double candY = p.y + dy;
140 if (!isInsideBounds(candX, candY, costmap))
142 if (!withinDeviation(candX, candY))
144 double candDist = queryDistance(wrapper, candX, candY);
145 if (candDist > bestDist)
160 if (bestDist >= clearance)
166 if (queryDistance(wrapper, p.x, p.y) <= 0.0)
168 const auto vertex = costmap.findClosestCollisionFreeVertex(
169 Eigen::Vector2f{
static_cast<float>(origX),
static_cast<float>(origY)},
170 static_cast<float>(maxDeviation));
172 if (vertex && withinDeviation(vertex->position.x(), vertex->position.y()))
174 p.x = vertex->position.x();
175 p.y = vertex->position.y();
181 bool validateTrajectory(
const std::vector<CenterPoint>& traj,
const std::string& label)
184 for (
int i = 0; i < static_cast<int>(traj.size()); ++i)
190 <<
"] : x=" << traj[i].x <<
" y=" << traj[i].y
191 <<
" theta=" << traj[i].theta <<
" v=" << traj[i].v;
198 ceres::Solver::Options makeSolverOptions(
const Params& params,
int maxIterations)
200 ceres::Solver::Options options;
201 options.max_num_iterations = maxIterations;
202 options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
203 options.minimizer_progress_to_stdout =
false;
204 options.num_threads = params.num_threads;
205 options.function_tolerance = 1e-6;
206 options.parameter_tolerance = 1e-8;
207 options.trust_region_strategy_type = ceres::LEVENBERG_MARQUARDT;
208 options.initial_trust_region_radius = params.initial_trust_region_radius;
209 options.max_trust_region_radius = params.max_trust_region_radius;
210 options.min_trust_region_radius = params.min_trust_region_radius;
211 options.min_relative_decrease = 1e-4;
212 options.use_nonmonotonic_steps =
true;
219 double w_pose_smooth;
222 double w_safe_tracking;
232 wrapper(&x, &y, &d_raw);
236 double obs_dist = d_raw - clearance;
237 if (obs_dist > max_dist)
241 double penetration = -obs_dist;
242 double max_pen = 3.0 * max_dist;
243 if (penetration > max_pen)
244 penetration = max_pen;
245 double ratio = penetration / max_dist;
246 return std::exp(2.0 * ratio) - 1.0 + 0.1 * ratio;
248 return (max_dist - obs_dist) / max_dist;
251 void logResidualBreakdown(
const std::vector<CenterPoint>& traj,
253 const std::vector<CenterPoint>& originalTargets,
254 const std::vector<CenterPoint>& safeTargets,
255 const Params& params,
256 const std::string& label)
258 const int N =
static_cast<int>(traj.size());
262 double c_obs_wp = 0.0;
263 double c_obs_seg = 0.0;
264 double c_smooth = 0.0;
265 double c_tracking = 0.0;
266 double c_safe_tracking = 0.0;
267 double c_spacing = 0.0;
270 for (
int i = 0; i < N; ++i)
272 double r = evalObstacleCost(wrapper, traj[i].x, traj[i].y, params.obs_max_distance, params.clearance);
277 for (
int i = 0; i < N - 1; ++i)
279 for (
int s = 1;
s <= params.segment_obstacle_samples; ++
s)
281 double t =
static_cast<double>(
s) / (params.segment_obstacle_samples + 1);
282 double mx = (1.0 - t) * traj[i].x + t * traj[i + 1].x;
283 double my = (1.0 - t) * traj[i].y + t * traj[i + 1].y;
284 double r = evalObstacleCost(wrapper, mx, my, params.obs_max_distance, params.clearance);
290 for (
int i = 1; i <= N - 2; ++i)
292 double rx = traj[i - 1].x - 2.0 * traj[i].x + traj[i + 1].x;
293 double ry = traj[i - 1].y - 2.0 * traj[i].y + traj[i + 1].y;
294 c_smooth += rx * rx + ry * ry;
298 for (
int i = 1; i <= N - 2; ++i)
300 double dx = traj[i].x - safeTargets[i].x;
301 double dy = traj[i].y - safeTargets[i].y;
302 double d = std::hypot(dx, dy);
303 double excess = d - params.tracking_deadzone;
306 if (excess <= params.tracking_max_deviation - params.tracking_deadzone)
308 double r = params.w_tracking_in_bounds * excess;
313 double over = excess - (params.tracking_max_deviation - params.tracking_deadzone);
314 double r = params.w_tracking_in_bounds * (params.tracking_max_deviation - params.tracking_deadzone)
315 + params.w_tracking_hard * over;
322 for (
int i = 1; i <= N - 2; ++i)
324 double rx = traj[i].x - originalTargets[i].x;
325 double ry = traj[i].y - originalTargets[i].y;
327 double w = params.pass3_w_tracking;
328 c_safe_tracking += (w * rx) * (w * rx) + (w * ry) * (w * ry);
332 double d_avg = computeAverageSpacing(traj);
335 for (
int i = 0; i < N - 1; ++i)
337 double dx = traj[i + 1].x - traj[i].x;
338 double dy = traj[i + 1].y - traj[i].y;
339 double dist_sq = dx * dx + dy * dy;
340 double r = params.pass3_w_spacing * (dist_sq - d_avg * d_avg);
344 double total = c_obs_wp + c_obs_seg + c_smooth + c_tracking + c_safe_tracking + c_spacing;
346 ARMARX_INFO <<
"=== Residual breakdown [" << label <<
"] ===";
347 ARMARX_INFO <<
"Obstacle (waypoints): " << std::scientific << std::setprecision(3) << c_obs_wp;
348 ARMARX_INFO <<
"Obstacle (segments): " << std::scientific << std::setprecision(3) << c_obs_seg;
349 ARMARX_INFO <<
"Position smoothness: " << std::scientific << std::setprecision(3) << c_smooth;
350 ARMARX_INFO <<
"Bounded tracking: " << std::scientific << std::setprecision(3) << c_tracking;
351 ARMARX_INFO <<
"Soft tracking: " << std::scientific << std::setprecision(3) << c_safe_tracking;
352 ARMARX_INFO <<
"Spacing: " << std::scientific << std::setprecision(3) << c_spacing;
353 ARMARX_INFO <<
"TOTAL: " << std::scientific << std::setprecision(3) << total;
356 void runPositionPass(std::vector<CenterPoint>& traj,
358 const std::vector<CenterPoint>& originalTargets,
359 const std::vector<CenterPoint>& safeTargets,
360 const PassWeights& w,
361 const Params& params,
364 const int N =
static_cast<int>(traj.size());
367 ARMARX_INFO <<
"[SPFA smoothing] need at least 4 points for position pass";
371 normalizeTrajectoryAngles(traj);
374 for (
int i = 1; i < N - 1; ++i)
376 double dx = traj[i].x - traj[i - 1].x;
377 double dy = traj[i].y - traj[i - 1].y;
378 if (std::hypot(dx, dy) < 1.0)
380 double ndx = (i + 1 < N) ? (traj[i + 1].x - traj[i].x) : 1.0;
381 double ndy = (i + 1 < N) ? (traj[i + 1].y - traj[i].y) : 0.0;
382 double nlen = std::hypot(ndx, ndy);
389 traj[i].x += 1.0 * ndx / nlen;
390 traj[i].y += 1.0 * ndy / nlen;
394 ceres::Problem problem;
395 for (
int i = 0; i < N; ++i)
396 problem.AddParameterBlock(&traj[i].x, 4);
398 problem.SetParameterBlockConstant(&traj.front().x);
399 problem.SetParameterBlockConstant(&traj.back().x);
402 for (
int i = 1; i < N - 1; ++i)
404 std::vector<int> constant_indices = {2, 3};
409 for (
int i = 1; i <= N - 2; ++i)
411 problem.AddResidualBlock(
412 new ceres::AutoDiffCostFunction<PositionSmoothResidual, 2, 4, 4, 4>(
421 double d_avg = computeAverageSpacing(traj);
424 for (
int i = 0; i < N - 1; ++i)
426 problem.AddResidualBlock(
427 new ceres::AutoDiffCostFunction<SpacingResidual, 1, 4, 4>(
435 for (
int i = 0; i < N; ++i)
437 problem.AddResidualBlock(
438 new ceres::AutoDiffCostFunction<ObstacleResidual2D, 1, 4>(
445 for (
int i = 0; i < N - 1; ++i)
447 for (
int s = 1;
s <= params.segment_obstacle_samples; ++
s)
449 double t =
static_cast<double>(
s) / (params.segment_obstacle_samples + 1);
450 problem.AddResidualBlock(
451 new ceres::AutoDiffCostFunction<SegmentObstacleResidual2D, 1, 4, 4>(
460 for (
int i = 1; i <= N - 2; ++i)
462 problem.AddResidualBlock(
463 new ceres::AutoDiffCostFunction<BoundedTrackingResidual, 1, 4>(
466 params.tracking_deadzone,
467 params.tracking_max_deviation,
468 params.w_tracking_in_bounds,
469 params.w_tracking_hard)),
475 if (w.w_tracking > 0.0)
477 for (
int i = 1; i <= N - 2; ++i)
479 problem.AddResidualBlock(
480 new ceres::AutoDiffCostFunction<SoftTrackingResidual, 2, 4>(
487 ceres::Solver::Options options = makeSolverOptions(params, maxIterations);
488 ceres::Solver::Summary summary;
489 ceres::Solve(options, &problem, &summary);
491 normalizeTrajectoryAngles(traj);
494 ARMARX_INFO <<
"[Position pass] Final cost: " << summary.final_cost
495 <<
" iterations: " << summary.num_successful_steps
496 <<
" termination: " << summary.termination_type;
499 bool repairWaypoints(std::vector<CenterPoint>& traj,
500 const algorithms::Costmap& costmap,
502 const std::vector<CenterPoint>& originalTargets,
503 const Params& params)
505 const int N =
static_cast<int>(traj.size());
507 for (
int i = 1; i < N - 1; ++i)
509 double d = queryDistance(wrapper, traj[i].x, traj[i].y);
510 if (d > params.clearance)
513 const auto vertex = costmap.findClosestCollisionFreeVertex(
514 Eigen::Vector2f{
static_cast<float>(originalTargets[i].x),
515 static_cast<float>(originalTargets[i].y)},
516 static_cast<float>(params.tracking_max_deviation));
521 const double dx = vertex->position.x() - originalTargets[i].x;
522 const double dy = vertex->position.y() - originalTargets[i].y;
523 if (std::hypot(dx, dy) > params.tracking_max_deviation + 1e-3)
526 traj[i].x = vertex->position.x();
527 traj[i].y = vertex->position.y();
533 double computeMaxDeviation(
const std::vector<CenterPoint>& traj,
534 const std::vector<CenterPoint>& originalTargets)
537 for (std::size_t i = 0; i < traj.size(); ++i)
539 double dx = traj[i].x - originalTargets[i].x;
540 double dy = traj[i].y - originalTargets[i].y;
541 maxDev = std::max(maxDev, std::hypot(dx, dy));
550 trajectory(trajectory),
560 std::vector<CenterPoint> traj = toCenterTrajectory(trajectory);
561 const std::vector<CenterPoint> originalTargets = traj;
563 if (!validateTrajectory(traj,
"pre-hill-climb"))
565 ARMARX_WARNING <<
"[SPFA smoothing] Invalid trajectory before pre-processing. Aborting.";
566 return {std::nullopt, std::nullopt,
false};
569 pushWaypointsToSafeZone(traj,
574 params.hill_climb_step_size,
575 params.hill_climb_max_iterations,
576 params.tracking_max_deviation);
578 const std::vector<CenterPoint> safeTargets = traj;
580 if (!validateTrajectory(traj,
"post-hill-climb"))
582 ARMARX_WARNING <<
"[SPFA smoothing] Invalid trajectory after pre-processing. Aborting.";
583 return {std::nullopt, std::nullopt,
false};
586 logResidualBreakdown(traj, wrapper, originalTargets, safeTargets, params,
"after pre-processing");
589 const std::array<PassWeights, 3> passes = {
590 PassWeights{params.pass1_w_obs,
591 params.pass1_w_pose_smooth,
592 params.pass1_w_spacing,
593 params.pass1_w_tracking,
594 params.pass1_w_safe_tracking},
595 PassWeights{params.pass2_w_obs,
596 params.pass2_w_pose_smooth,
597 params.pass2_w_spacing,
598 params.pass2_w_tracking,
599 params.pass2_w_safe_tracking},
600 PassWeights{params.pass3_w_obs,
601 params.pass3_w_pose_smooth,
602 params.pass3_w_spacing,
603 params.pass3_w_tracking,
604 params.pass3_w_safe_tracking}};
606 for (
size_t pass = 0; pass < passes.size(); ++pass)
608 ARMARX_INFO <<
"[SPFA smoothing] Starting position pass " << (pass + 1);
609 runPositionPass(traj, wrapper, originalTargets, safeTargets, passes[pass], params, params.max_iterations);
610 logResidualBreakdown(traj, wrapper, originalTargets, safeTargets, params,
611 "after position pass " + std::to_string(pass + 1));
613 if (!validateTrajectory(traj,
"post-position-pass-" + std::to_string(pass + 1)))
615 ARMARX_WARNING <<
"[SPFA smoothing] NaN/Inf after position pass " << (pass + 1)
616 <<
". Reverting to pre-processed safe path.";
622 normalizeTrajectoryAngles(traj);
627 bool collisionFree = checker.
check(smoothedTrajectory);
631 ARMARX_WARNING <<
"[SPFA smoothing] Initial smoothing has collisions. Running repair pass.";
633 if (!repairWaypoints(traj, costmap, wrapper, originalTargets, params))
635 ARMARX_ERROR <<
"[SPFA smoothing] Repair failed: could not project colliding waypoints "
636 "to collision-free positions within max deviation.";
640 PassWeights repairWeights{params.repair_w_obs,
641 params.repair_w_pose_smooth,
642 params.repair_w_spacing,
643 params.repair_w_tracking,
644 params.repair_w_safe_tracking};
645 runPositionPass(traj, wrapper, originalTargets, safeTargets, repairWeights, params, params.repair_max_iterations);
646 logResidualBreakdown(traj, wrapper, originalTargets, safeTargets, params,
"after repair pass");
648 smoothedTrajectory = toGlobalTrajectory(traj);
649 collisionFree = checker.
check(smoothedTrajectory);
653 const double maxDeviation = computeMaxDeviation(traj, originalTargets);
654 const bool deviationOk = maxDeviation <= params.tracking_max_deviation + 1e-3;
656 if (collisionFree && deviationOk)
658 ARMARX_INFO <<
"[SPFA smoothing] Smoothed trajectory is collision-free. Max deviation: "
659 << maxDeviation <<
" mm";
664 ARMARX_WARNING <<
"[SPFA smoothing] Smoothed trajectory still has collisions after repair.";
666 ARMARX_WARNING <<
"[SPFA smoothing] Max deviation too large: " << maxDeviation
667 <<
" mm (limit: " << params.tracking_max_deviation <<
" mm).";
670 for (
int i = 0; i < static_cast<int>(traj.size()); ++i)
672 ARMARX_DEBUG << std::fixed << std::setprecision(3) <<
"i=" << i <<
" : x=" << traj[i].x
673 <<
" y=" << traj[i].y <<
" theta=" << traj[i].theta <<
" v=" << traj[i].v;
677 const bool success = collisionFree && deviationOk;
678 return {smoothedTrajectory, preprocessedTrajectory,
success};
Differentiable wrapper around the standard 2-D distance-to-obstacle costmap.
OptimizationResult optimize()
SPFASmoothing(const core::GlobalTrajectory &trajectory, const algorithms::Costmap &costmap, const Params ¶ms)
io::SmoothingParams Params
bool check(const core::GlobalTrajectory &trajectory, bool logDetails=false) const
#define ARMARX_INFO
The normal logging level.
#define ARMARX_ERROR
The logging level for unexpected behaviour, that must be fixed.
#define ARMARX_DEBUG
The logging level for output that is only interesting while debugging.
#define ARMARX_WARNING
The logging level for unexpected behaviour, but not a serious problem.
double s(double t, double s0, double v0, double a0, double j)
std::vector< CenterPoint > toCenterTrajectory(const armarx::navigation::core::GlobalTrajectory >raj)
armarx::navigation::core::GlobalTrajectory toGlobalTrajectory(const std::vector< CenterPoint > &traj)
This file is part of ArmarX.
std::multimap< std::string, std::string > Params
bool isfinite(const std::vector< T, Ts... > &v)
ceres::SubsetParameterization SubsetConstraint
This file is part of ArmarX.
#define SET_CONSTRAINT(problem, param, constraint)
Bounded tracking residual.