轨迹优化 | 图解欧氏距离场与梯度场算法(附ROS C++/Python实现)

2024-06-17 12:20

本文主要是介绍轨迹优化 | 图解欧氏距离场与梯度场算法(附ROS C++/Python实现),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

目录

  • 0 专栏介绍
  • 1 什么是距离场?
  • 2 欧氏距离场计算原理
  • 3 双线性插值与欧式梯度场
  • 4 仿真实现
    • 4.1 ROS C++实现
    • 4.2 Python实现

0 专栏介绍

🔥课程设计、毕业设计、创新竞赛、学术研究必备!本专栏涉及更高阶的运动规划算法实战:曲线生成与轨迹优化、碰撞模型与检测、多智能体群控、深度强化学习运动规划、社会性导航、全覆盖路径规划等内容,每个模型都包含代码实现加深理解。

🚀详情:运动规划实战进阶


1 什么是距离场?

距离场(Distance Field)也称为距离变换(Distance Transform),在图像处理和模式识别中是一种重要工具,其核心思想是将图像中每个像素点的值表示为到最近的目标像素的距离度量,更高维的情况依次类推。在计算机视觉领域,距离场是图像分割和配准的基础,在运动规划中则常用于地图建模(广义Voronoi图计算)和轨迹梯度优化。

在这里插入图片描述

给定点集 G \mathcal{G} G和其上的目标点集 Q ⊆ G Q\subseteq \mathcal{G} QG,可以定义函数

D I ( p ) = min ⁡ q ∈ G ( d ( p , q ) + I ( q ) ) \mathcal{D} _I\left( p \right) =\min _{q\in \mathcal{G}}\left( d\left( p,q \right) +I\left( q \right) \right) DI(p)=qGmin(d(p,q)+I(q))

计算任意一点 p p p Q Q Q的最近距离,其中距离度量 d ( ⋅ , ⋅ ) d\left( \cdot ,\cdot \right) d(,)定义了距离场的属性,指示函数

I Q ( q ) = { 0 , i f q ∈ Q ∞ , o t h e r w i s e I_Q\left( q \right) =\begin{cases} 0, \mathrm{if} q\in Q\\ \infty , \mathrm{otherwise}\\\end{cases} IQ(q)={0,ifqQ,otherwise

常用的 d ( ⋅ , ⋅ ) d\left( \cdot ,\cdot \right) d(,)

  • 欧氏距离,此时称为欧氏距离变换(Euclidean Distance Transform, EDT)
  • 曼哈顿距离,此时称为曼哈顿距离变换(Manhattan Distance Transform, MDT)
  • 切比雪夫距离,此时称为切比雪夫距离变换(Chebyshev Distance Transform, CDT)

2 欧氏距离场计算原理

n n n维距离场可以通过一维距离场迭代计算得到,因此只需要讨论一维EDT的计算即可。如下图左侧所示为初始计算轴的EDT计算,右侧所示为更一般的情况,此时障碍物处的采样函数叠加了前轴计算信息。

在这里插入图片描述

注意到 O \mathcal{O} O定义了一系列以障碍物 q ∈ O q\in \mathcal{O} qO为顶点的抛物线,而 ∀ p ∈ G \forall p\in \mathcal{G} pG在系列抛物线形成的下包络的投影组成了距离场 D f ( p ) \mathcal{D} _f\left( p \right) Df(p),而下包络的计算与抛物线交点有关。联立两条抛物线 ( s − q ) 2 + f ( q ) = ( s − r ) 2 + f ( r ) \left( s-q \right) ^2+f\left( q \right) =\left( s-r \right) ^2+f\left( r \right) (sq)2+f(q)=(sr)2+f(r)可得

s = ( f ( r ) + r 2 ) − ( f ( q ) + q 2 ) 2 r − 2 q s=\frac{\left( f\left( r \right) +r^2 \right) -\left( f\left( q \right) +q^2 \right)}{2r-2q} s=2r2q(f(r)+r2)(f(q)+q2)

即任意两条抛物线有且仅有一个交点 s s s。设 K \mathcal{K} K为实际组成下包络的抛物线集合, v ( k ) v\left( k \right) v(k)表示其中第 k k k条抛物线的顶点, z ( k ) z\left( k \right) z(k)表示第 k k k条和第 k − 1 k-1 k1条抛物线的交点,区间 [ z ( k ) , z ( k + 1 ) ) \left[ z\left( k \right) ,z\left( k+1 \right) \right) [z(k),z(k+1))表示第 k k k条抛物线的下包络范围。在遍历求解下包络过程中,对于新的抛物线 e e e,其与 K \mathcal{K} K中最新的一条抛物线 k k k的交点 s s s有两种情况:

  • s > z ( k ) s>z\left( k \right) s>z(k),则将 e e e添加到 K \mathcal{K} K中并更新 v ( k ) v\left( k \right) v(k) z ( k ) z\left( k \right) z(k)
  • s ⩽ z ( k ) s\leqslant z\left( k \right) sz(k),则第 k k k条抛物线不参与构成下包络,应从 K \mathcal{K} K中删除并重新计算 e e e与新的第 k k k条抛物线的交点直至 s s s z ( k ) z(k) z(k)右侧;

算法流程如表所示

在这里插入图片描述

3 双线性插值与欧式梯度场

在求解梯度过程中,需要计算离散距离场的线性插值函数。以二维环境的双线性插值为例,设已知离散空间有四点坐标为 A ( x 1 , y 1 ) A\left( x_1,y_1 \right) A(x1,y1) B ( x 2 , y 1 ) B\left( x_2,y_1 \right) B(x2,y1) C ( x 1 , y 2 ) C\left( x_1,y_2 \right) C(x1,y2) D ( x 2 , y 2 ) D\left( x_2,y_2 \right) D(x2,y2),对其中任意一点 P P P,首先在 x x x方向上插值

{ f ( R 1 ) = f ( x , y 1 ) = x 2 − x x 2 − x 1 f ( A ) + x − x 1 x 2 − x 1 f ( B ) f ( R 2 ) = f ( x , y 2 ) = x 2 − x x 2 − x 1 f ( C ) + x − x 1 x 2 − x 1 f ( D ) \begin{cases} f\left( R_1 \right) =f\left( x,y_1 \right) =\frac{x_2-x}{x_2-x_1}f\left( A \right) +\frac{x-x_1}{x_2-x_1}f\left( B \right)\\ f\left( R_2 \right) =f\left( x,y_2 \right) =\frac{x_2-x}{x_2-x_1}f\left( C \right) +\frac{x-x_1}{x_2-x_1}f\left( D \right)\\\end{cases} {f(R1)=f(x,y1)=x2x1x2xf(A)+x2x1xx1f(B)f(R2)=f(x,y2)=x2x1x2xf(C)+x2x1xx1f(D)

再基于插值点 R 1 R_1 R1 R 2 R_2 R2进行 y y y方向的插值

f ( P ) = f ( x , y ) = y 2 − y y 2 − y 1 f ( R 1 ) + y − y 1 y 2 − y 1 f ( R 2 ) f\left( P \right) =f\left( x,y \right) =\frac{y_2-y}{y_2-y_1}f\left( R_1 \right) +\frac{y-y_1}{y_2-y_1}f\left( R_2 \right) f(P)=f(x,y)=y2y1y2yf(R1)+y2y1yy1f(R2)

展开可得矩阵形式

f ( x , y ) = 1 ( x 2 − x 1 ) ( y 2 − y 1 ) [ x 2 − x x − x 1 ] T [ f ( A ) f ( C ) f ( B ) f ( D ) ] [ y 2 − y y − y 1 ] f\left( x,y \right) =\frac{1}{\left( x_2-x_1 \right) \left( y_2-y_1 \right)}\left[ \begin{array}{c} x_2-x\\ x-x_1\\\end{array} \right] ^T\left[ \begin{matrix} f\left( A \right)& f\left( C \right)\\ f\left( B \right)& f\left( D \right)\\\end{matrix} \right] \left[ \begin{array}{c} y_2-y\\ y-y_1\\\end{array} \right] f(x,y)=(x2x1)(y2y1)1[x2xxx1]T[f(A)f(B)f(C)f(D)][y2yyy1]

在这里插入图片描述

对于离散栅格坐标而言, x 2 − x 1 = y 2 − y 1 = 1 x_2-x_1=y_2-y_1=1 x2x1=y2y1=1,设 Δ x = x − x 1 \varDelta x=x-x_1 Δx=xx1 Δ y = y − y 1 \varDelta y=y-y_1 Δy=yy1,则插值函数简化为

f ( Δ x , Δ y ) = [ 1 − Δ x Δ x ] T [ f ( A ) f ( C ) f ( B ) f ( D ) ] [ 1 − Δ y Δ y ] f\left( \varDelta x,\varDelta y \right) =\left[ \begin{array}{c} 1-\varDelta x\\ \varDelta x\\\end{array} \right] ^T\left[ \begin{matrix} f\left( A \right)& f\left( C \right)\\ f\left( B \right)& f\left( D \right)\\\end{matrix} \right] \left[ \begin{array}{c} 1-\varDelta y\\ \varDelta y\\\end{array} \right] f(Δx,Δy)=[1ΔxΔx]T[f(A)f(B)f(C)f(D)][1ΔyΔy]

则函数 f f f在任意一点的梯度为

∇ f ( Δ x , Δ y ) = [ ∂ f ( Δ x , Δ y ) ∂ Δ x ∂ f ( Δ x , Δ y ) ∂ Δ y ] = [ [ − 1 1 ] T [ f ( A ) f ( C ) f ( B ) f ( D ) ] [ 1 − Δ y Δ y ] [ 1 − Δ x Δ x ] T [ f ( A ) f ( C ) f ( B ) f ( D ) ] [ − 1 1 ] ] \nabla f\left( \varDelta x,\varDelta y \right) =\left[ \begin{array}{c} \frac{\partial f\left( \varDelta x,\varDelta y \right)}{\partial \varDelta x}\\ \frac{\partial f\left( \varDelta x,\varDelta y \right)}{\partial \varDelta y}\\\end{array} \right] =\left[ \begin{array}{c} \left[ \begin{array}{c} -1\\ 1\\\end{array} \right] ^T\left[ \begin{matrix} f\left( A \right)& f\left( C \right)\\ f\left( B \right)& f\left( D \right)\\\end{matrix} \right] \left[ \begin{array}{c} 1-\varDelta y\\ \varDelta y\\\end{array} \right]\\ \left[ \begin{array}{c} 1-\varDelta x\\ \varDelta x\\\end{array} \right] ^T\left[ \begin{matrix} f\left( A \right)& f\left( C \right)\\ f\left( B \right)& f\left( D \right)\\\end{matrix} \right] \left[ \begin{array}{c} -1\\ 1\\\end{array} \right]\\\end{array} \right] f(Δx,Δy)=[Δxf(Δx,Δy)Δyf(Δx,Δy)]=[11]T[f(A)f(B)f(C)f(D)][1ΔyΔy][1ΔxΔx]T[f(A)f(B)f(C)f(D)][11]

4 仿真实现

4.1 ROS C++实现

核心代码如下所示

GradientLayer::updateCosts(nav2_costmap_2d::Costmap2D & master_grid, int min_i, int min_j,int max_i,int max_j)
{if (!enabled_) {return;}unsigned char * master_array = master_grid.getCharMap();unsigned int size_x = master_grid.getSizeInCellsX(), size_y = master_grid.getSizeInCellsY();// Fixing window coordinates with map size if necessary.min_i = std::max(0, min_i);min_j = std::max(0, min_j);max_i = std::min(static_cast<int>(size_x), max_i);max_j = std::min(static_cast<int>(size_y), max_j);// Simply computing one-by-one cost per each cellint gradient_index;for (int j = min_j; j < max_j; j++) {// Reset gradient_index each time when reaching the end of re-calculated window// by OY axis.gradient_index = 0;for (int i = min_i; i < max_i; i++) {int index = master_grid.getIndex(i, j);// setting the gradient costunsigned char cost = (LETHAL_OBSTACLE - gradient_index*GRADIENT_FACTOR)%255;if (gradient_index <= GRADIENT_SIZE) {gradient_index++;} else {gradient_index = 0;}master_array[index] = cost;}}
}

4.2 Python实现

欧氏距离场核心代码

def compute(self, f_get, mat: np.ndarray, dim: int) -> np.array:"""Compute distance field along one-dimension baseon sample funciton f_get.Parameters:f_get (function): sample functionmat (np.ndarray): the matrix to transformdim (int): the dimension to transformReturns:df (np.array): the distance field along dim-dimension"""# initialzationk = 0n = mat.shape[dim]v, z = [0 for _ in range(n)], [0 for _ in range(n + 1)]z[0], z[1] = -self.INF, self.INF# envelopefor q in range(1, n):s = ((f_get(v[k]) + v[k] ** 2) - (f_get(q) + q ** 2)) / (2 * (v[k] - q))while s <= z[k]:k -= 1s = ((f_get(v[k]) + v[k] ** 2) - (f_get(q) + q ** 2)) / (2 * (v[k] - q))k += 1v[k] = qz[k], z[k + 1] = s, self.INF# distance calculationk = 0edf = np.zeros((n, ))for q in range(n):while z[k + 1] < q:k += 1edf[q] = (q - v[k]) ** 2 + f_get(v[k])return edf

欧氏梯度场核心代码:

def gradient(self, df: np.ndarray, x: float, y: float) -> np.array:"""To obtain the gradient at (x, y) in the distance field through bilinear interpolation.Parameters:df (np.ndarray): the distance fieldx/y (float): the query coordinateReturns:g(x, y): the gradient at (x, y)"""m, n = df.shapex, y = max(min(n - 1, x), 0), max(min(m - 1, y), 0)xi, yi = int(x), int(y)dx, dy = x - xi, y - yixi, yi = max(min(n - 1, xi), 0), max(min(m - 1, yi), 0)xp, yp = max(min(n - 1, xi + 1), 0), max(min(m - 1, yi - 1), 0)bl, br = df[yi, xi], df[yi, xp]tl, tr = df[yp, xi], df[yp, xp]return np.array([(1 - dy) * (br - bl) + dy * (tr + tl),-((1 - dx) * (tl - bl) + dx * (tr - br))])

效果如下所示
在这里插入图片描述

一个实际地图的案例如下

在这里插入图片描述

完整工程代码请联系下方博主名片获取


🔥 更多精彩专栏

  • 《ROS从入门到精通》
  • 《Pytorch深度学习实战》
  • 《机器学习强基计划》
  • 《运动规划实战精讲》

👇源码获取 · 技术交流 · 抱团学习 · 咨询分享 请联系👇

这篇关于轨迹优化 | 图解欧氏距离场与梯度场算法(附ROS C++/Python实现)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1069424

相关文章

pandas中位数填充空值的实现示例

《pandas中位数填充空值的实现示例》中位数填充是一种简单而有效的方法,用于填充数据集中缺失的值,本文就来介绍一下pandas中位数填充空值的实现,具有一定的参考价值,感兴趣的可以了解一下... 目录什么是中位数填充?为什么选择中位数填充?示例数据结果分析完整代码总结在数据分析和机器学习过程中,处理缺失数

Golang HashMap实现原理解析

《GolangHashMap实现原理解析》HashMap是一种基于哈希表实现的键值对存储结构,它通过哈希函数将键映射到数组的索引位置,支持高效的插入、查找和删除操作,:本文主要介绍GolangH... 目录HashMap是一种基于哈希表实现的键值对存储结构,它通过哈希函数将键映射到数组的索引位置,支持

Pandas使用AdaBoost进行分类的实现

《Pandas使用AdaBoost进行分类的实现》Pandas和AdaBoost分类算法,可以高效地进行数据预处理和分类任务,本文主要介绍了Pandas使用AdaBoost进行分类的实现,具有一定的参... 目录什么是 AdaBoost?使用 AdaBoost 的步骤安装必要的库步骤一:数据准备步骤二:模型

使用Pandas进行均值填充的实现

《使用Pandas进行均值填充的实现》缺失数据(NaN值)是一个常见的问题,我们可以通过多种方法来处理缺失数据,其中一种常用的方法是均值填充,本文主要介绍了使用Pandas进行均值填充的实现,感兴趣的... 目录什么是均值填充?为什么选择均值填充?均值填充的步骤实际代码示例总结在数据分析和处理过程中,缺失数

Java对象转换的实现方式汇总

《Java对象转换的实现方式汇总》:本文主要介绍Java对象转换的多种实现方式,本文通过实例代码给大家介绍的非常详细,对大家的学习或工作具有一定的参考借鉴价值,需要的朋友参考下吧... 目录Java对象转换的多种实现方式1. 手动映射(Manual Mapping)2. Builder模式3. 工具类辅助映

如何使用 Python 读取 Excel 数据

《如何使用Python读取Excel数据》:本文主要介绍使用Python读取Excel数据的详细教程,通过pandas和openpyxl,你可以轻松读取Excel文件,并进行各种数据处理操... 目录使用 python 读取 Excel 数据的详细教程1. 安装必要的依赖2. 读取 Excel 文件3. 读

Go语言开发实现查询IP信息的MCP服务器

《Go语言开发实现查询IP信息的MCP服务器》随着MCP的快速普及和广泛应用,MCP服务器也层出不穷,本文将详细介绍如何在Go语言中使用go-mcp库来开发一个查询IP信息的MCP... 目录前言mcp-ip-geo 服务器目录结构说明查询 IP 信息功能实现工具实现工具管理查询单个 IP 信息工具的实现服

Python的time模块一些常用功能(各种与时间相关的函数)

《Python的time模块一些常用功能(各种与时间相关的函数)》Python的time模块提供了各种与时间相关的函数,包括获取当前时间、处理时间间隔、执行时间测量等,:本文主要介绍Python的... 目录1. 获取当前时间2. 时间格式化3. 延时执行4. 时间戳运算5. 计算代码执行时间6. 转换为指

利用Python调试串口的示例代码

《利用Python调试串口的示例代码》在嵌入式开发、物联网设备调试过程中,串口通信是最基础的调试手段本文将带你用Python+ttkbootstrap打造一款高颜值、多功能的串口调试助手,需要的可以了... 目录概述:为什么需要专业的串口调试工具项目架构设计1.1 技术栈选型1.2 关键类说明1.3 线程模

SpringBoot基于配置实现短信服务策略的动态切换

《SpringBoot基于配置实现短信服务策略的动态切换》这篇文章主要为大家详细介绍了SpringBoot在接入多个短信服务商(如阿里云、腾讯云、华为云)后,如何根据配置或环境切换使用不同的服务商,需... 目录目标功能示例配置(application.yml)配置类绑定短信发送策略接口示例:阿里云 & 腾