车辙轨迹预测线算法讲解和代码实现
附赠自动驾驶学习资料和量产经验:链接
一、效果展示
鱼眼
去畸变
畸变矫正,正视图
畸变矫正,正视图,更广的角度;还是鱼眼图看得多一点
看下市场上的:
实验看起来算一致。
视频:
00:08
00:35
00:36
00:32
下面介绍原理。
二、原理介绍
两部分:分为车辙点生成和点的视图转换
1.车辙点生成
条件:(1)在AVM鸟瞰图上;(2)低速;
工具:(1)横向运动学模型;(2)在(1)基础上衍生出的小技巧,更简单
额外:(1)刻度线,1米2米3米;(2)车体,轮子,各种线
调优:(1)车辙点稀疏化;(2)用小技巧调整计算方式
2.点的视图间转换
视图类型:(1)鸟瞰图BEV,也即AVM;(2)鱼眼图FEV;(3)去畸变图(鱼眼图去畸变)(4)去畸变扳正图,把去畸变中前方物体扳正,正视图,相当于摄像头垂直拍车子前面和垂直拍车子后面;(5)更大视野范围的去畸变图
这些不同视图,统统以(1)鸟瞰图下的车辙点进行转换。
条件:(1)鸟瞰图上车辙点
调优:(1)这些视图上的车辙点稀疏化
分一条龙服务,从头到尾所有涉及一遍,再之后,直接调用现成的AVM标定好以后的结果计算车辙。所以,先从复杂的开始讲,再到简单的方式;
2.1 车辙点生成
2.1.1 预备知识:旋转平移
点旋转 和 坐标系旋转
点旋转
坐标系 绕坐标原点 旋转(逆时针)
坐标系旋转
2.1.2 Lateral Kinematics Model 横向运动学模型
Lateral Vehicle Dynamics
在书的第21页,pdf的第48页
在书的第22、23、24页,pdf的49、50、51
这个模型的关键点在于:
神奇的假设
既然假设了 轨道运动 ,那么是圆周运动,这个是瞬时的,如果我们迭代计算,下一时刻是否还符合这个运动呢? 难以理解,只能实车看一下吧,反正是低速,我们需要在视频里打点,看低速恒方向盘转角,的轨迹,N多个瞬时时刻(dt)后是否能到那里。我没有试,但从录得视频来看,差不多。我们假定就是这样的理想轨迹,此时,就是转圈圈儿,预测轨迹就是转圈圈儿,好像很扯。
恒速恒角理想预测
这里需要说明的是,dt很小,即每次迭代间的delta_time取得非常小,如取1ms,3ms,大了的话,这个模型要失效。
dt (delta_time) 不宜过大,当然也有技巧方法,即小dt迭代次数多几次再画点
我们选取计算点在A点上,以车身方向为X轴,是一个相对坐标系,模型简化为
若取 O-XY,水平方向为X轴,垂直方向为Y轴
最后,以车机实际显示,建立坐标系,车辆始终被保持在车机显示的正中央,车身方向就是Y轴,垂直方向
这个模型大有用处,在AVM中,不仅用来生成车辙,也可用来做透明底盘
其实看到这里,就大概明朗了,用这个 ↑ 做迭代算就行了。
一言以蔽之:
更好用的是,直接用δ,计算轨道圆的方程,因为在车上的车机显示板上,车始终在屏幕正中央,这很关键。直接计算一个圆的方程,大大化简所有计算。
理想情况下,就是在一个圆上做圆周运动。
// (X - X0)^2 + (Y - Y0)^2 = R^2
// X^2 - 2*X0*X + X0^2 + (Y - Y0)^2 - R^2 = 0
// a*x^2 + b*x + c = 0
// a = 1, b = -2*x0, c = ..
// delta = b^2 - 4ac
// X = , Y = ;一元二次方程的求根公式:x=[-b±√(b²-4ac)]/2a
// X0,Y0就是圆心,同心圆
2.1.3 点生成
复杂的方法:
pipeline
-
从运动学模型,计算出 选取点 的坐标平移dx,dy,和车身的角度偏转dθ,
-
从点旋转Rt和{dx,dy,dθ},计算出t(平移向量)
-
计算坐标轴旋转平移的R’t’,可以发现,t’实际等于t
-
通过上述算得的条件,迭代计算每一次更新后的坐标
需要注意的是,坐标系的选取的原点,以车中心为坐标系原点,
此原点既拿来(1)经过点旋转;和也拿来(2)经过坐标系旋转。一个点用两次。作为运动学模型和坐标系旋转之间的桥梁。
伪代码在下章节。
看代码吧,复杂的方法太蹩脚了,不值得讲。倒是点的视图间转换复杂的方法值得讲。
简单的方法:
简单的方法就是
极致的方法:
用一元二次求根公式即可,
double c = std::pow(circle_center_x, 2) +
std::pow(y - circle_center_y, 2) - std::pow(_R, 2);
double delta = std::pow(b, 2) - 4. * a * c;
float x;
if (delta >= 0) {
int coef = direction == 0 ? -1 : 1;
x = float(-b + coef * std::sqrt(delta)) / 2.f;
} else {
continue;
}
还是代码说好说一点。
x根肯定要为正,负的不要。
2.2 点的视图间转换
两种方式:(1)费时费力的处理:通过remap warp;(2)简易的换算:直接算
鸟瞰图,到 去畸变图,存在一个Homography,就是perspective transform,透视变换
去畸变图,到 FEV鱼眼图,存在一个 去畸变模型的查找表。(我在这里卡了很久)
原因在于,去畸变图 反映射回 畸变图,这个 这个本身就不对的道理,我们不要拿图去反映射,
而是拿points of rut trajectory prediction车辙轨迹预测的点的坐标,通过一吨操作,在鱼眼图上去画,并不是拿画好后的BEV(AVM鸟瞰图)去做变换,想通这里,接下来就好办了。
加粗一下:
得到鸟瞰图到去畸变图的关系后,无论是接下来到:(a)去畸变扳正 正视图、(b)去畸变扳正更大视野正视图、(c)鱼眼图、(d)其他任意按需设计的视图, 都只需从去畸变图出发,去找对应的Homograpy单应矩阵,然后把点从去畸变投到相应视图上。
单应矩阵是什么?:
看过的,参考下吧:
单应性矩阵的理解及求解_机器视觉001的博客-CSDN博客_单应性矩阵
利用逆透视变换获取车载图象的鸟瞰图_玉米味土豆片的博客-CSDN博客_opencv 鸟瞰图
水木十三:opencv 单应性矩阵计算函数findHomography
findHomography()函数详解_奔跑的小木的博客-CSDN博客_findhomography
getPerspectiveTransform()与findHomography()与 HomographyBasedEstimator()
warpPerspective函数_普通网友的博客-CSDN博客_warpperspective
Rodrigues的推导: 罗德里格斯公式Rodrigues’Rotation Formula推导
2.2.1 鱼眼去畸变模型
我们先看看 去畸变,(由畸变图到去畸变), OpenCV文档:
去畸变模型
def calDistCoord(coords_undist, dist_coeffs, w_new, h_new):
"""
对opencv去畸变模型的代码实现
输入:
coords_undist,去畸变的坐标索引,相当于(无畸变)入射光线的坐标
dist_coeffs,畸变参数
w_new, h_new,去畸变图的大小
"""
K = np.array([[fx,0,cx], [0,fy,cy], [0,0,1]]) # cx = w/2, cy = h/2
K_new = np.array([[fx,0,w_new/2], [0,fy,h_new/2], [0,0,1]])
# 用矩阵的方式来做
coords_undist[:, :, 0] = (coords_undist[:, :, 0] - K_new[0, 2]) / K_new[0, 0]
coords_undist[:, :, 1] = (coords_undist[:, :, 1] - K_new[1, 2]) / K_new[1, 1]
r_undist = np.sqrt(point_undist[:, :, 0] ** 2 + coords_undist[:, :, 1] ** 2)[
:, :, np.newaxis
] # 齐次坐标
angle_undist = np.arctan(r_undist) # / np.pi * 180
r_dist = func(angle_undist, *dis_coeffs) # opencv里是θd,这里是r_dist,不用opencv的表示
coords_dist = np.zeros_like(coords_undist)
index = r_undist[:, :, 0] != 0 # r_undist[cy]cx]等于0,坐标原点,这个拟射到这个点的光线没折射
# 通过矩阵的形式来做
coords_dist[index] = r_dist[index] / r_undist[index] * point_undist[index]
# 算完以后,还要 u = x'*fx + cx_new, v = y'*fy + cy_new
coords_dist[:, :, 0] = coords_dist[:, :, 0] * K[0, 0] + K[0, 2]
coords_dist[:, :, 1] = coords_dist[:, :, 1] * K[1, 1] + K[1, 2]
return coords_dist
scale = 2
point_undist = np.meshgrid(np.arange(640 * scale), np.arange(480 * scale))
point_undist = np.stack(point_undist, axis=2).astype(np.float32) # h, w, 2 存的坐标
img_size = (640 * scale, 480 * scale) # W, H
map = calDistCoord(point_undist, *img_size) # 这个是opencv的fisheye去畸变模型实现
// 相机内参已知,畸变参数已知
# map.shape = (960, 1280, 2) (h, w, xy)
# 举例
x, y = 640, 480 # 原点,鱼眼相机小孔成像的直射点,此方向光线无畸变
new_x = map[y, x, 0] # = 640 map_x(x,y) = map[:,:,0]
new_y = map[y, x, 1] # = 480 map_y(x,y) = map[:,:,1]
x, y = 640 - 640, 480 - 480
new_x = map[0, 0, 0] # = 299.0378
new_y = map[0, 0, 1] # = 224.27834
x, y = 639 + 640, 479 + 480
new_x = map[y, x, 0] # = 980.82605
new_y = map[y, x, 1] # = 735.4862
print(new_x, new_y) # 这个新坐标就是鱼眼图上的 车辙轨迹点
# 直接在map里面找
# 需要重点说明的是,这个↓是 鱼眼图去畸变,得到去畸变后的图,我们在remap前就停止操作
# undist_img = cv2.remap(dist_img, map[:, :, 0], map[:, :, 1], cv2.INTER_LINEAR)
# 不进行这一步 !↑
# # remap: dst(x,y)=src(map_x(x,y),map_y(x,y))
# # 题外话:
# # remap: Applies a generic geometrical transformation to an image.
# # remap真是个神奇的函数
""" cv2.remap(src_img, map1, map2, interpolation)
src_img:原始图像
map1:目标图像中的每一个像素都对应于(映射于)原始图像的中的某个像素,
将对应的原始图像的像素坐标的col存放在map1中。map1 的shape和目标图像一样。
map2:目标图像中的每一个像素都对应于(映射于)原始图像的中的某个像素,
将对应的原始图像的像素坐标的row存放在map2中。map2 的shape和目标图像一样。
interpolation: 插值方式,不支持 INTER_AREA
"""
# 明明就是 畸变模型,竟然通过remap去畸变了
# https://github.com/opencv/opencv/blob/4.6.0/modules/imgproc/src/imgwarp.cpp#L1664
# meshgrid我一直搞得很迷糊,这里不记录下,不占用点篇幅:
opencv这个fisheye里写的模型,
Fisheye distortion
那应该是:
Juho Kannala and Sami Brandt. A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses.IEEE transactions on pattern analysis and machine intelligence, 28:1335–40, 09 2006.
用的这篇文章
OpenCV这样写。
再缕下思路: 通过 标定布的点(bev鸟瞰图视角) 与去畸变图上对应的点,
算得 鸟瞰图 ->到 去畸变图的 H。 把之前 鸟瞰图上的 轨迹预测点 通过H 算到 去畸变图上,
再算到 畸变图(鱼眼图)上,整个流程就走完了。
那其实可以直接 算 从鸟瞰图到 鱼眼图的映射吗? 没有。用去畸变图中转一下。
2.2.2 去畸变图扳正成正视图
这里一开始有个误区:
如果用AVM的流程直接调用,拿鱼眼图整张图变换到 扳正后的去畸变图,这里会使得车辙点被插值,这样要不得
要不得
正确的做法
复杂的方法是 算从去畸变到 扳正后的正视图的 map, map_x 和 map_y, 然后,输入BEV鸟瞰图下的坐标,得到 正视图的坐标。
这里需要算两个Homograpy单应性矩阵: H_bev2undist, H_undist2thisview,
,然后,用两个Homo 去计算map, 把map存为bin文件,用的时候直接调。麻烦啊,初始化占时间
简单的方法是
用两个Homo,把输入坐标直接用Homo做透视变换。方便,
复杂的方法是走的 AVM的标定流程,计算车辙时不需要。而且甚至最终解决方案是:AVM标定时把 中间结果Homograpy给存下来。 这样 省时省力。
在下章代码实现里,也说一下,怎么生成map的bin文件。
2.2.3 更广的视角
无非是套娃,再加多点变换
2.2.4 其他
同上
三、代码实现
伪了
3.1 点生成
复杂的方法:
伪代码
float model_point[2] = {200, 200}; // 举例
float dth=0, dx=0, dy=0; // theta, x, y
rotate_center_init(...);
lateral_kinematics_model(model_point) // calculate: dth, dx, dy
void calc_trans_vector() {
// use lateral kinematics model calculation
//使用前轮模型,这里用前轮中心点坐标
Point3D object_point = m_front_wheel_center;
// 这里object_point的坐标需要以车中心为坐标系为原点,需要换算一下
// TODO ( )
// add the predictions from kinematics model
Point3D predicted_pt;
predicted_pt.x = object_point.x + delta_x;
predicted_pt.y = object_point.y + delta_y;
Point3D rotated_pt;
rotate(rotated_pt, object_point);
// get translation vector for wheel center
translation_vec[0] = predicted_pt.x - rotated_pt.x;
translation_vec[1] = predicted_pt.y - rotated_pt.y;
} // 这个是运用运动学模型计算车辙点的关键一招
float T[2] = {上↑面}; // 平移向量t'等于t
// R1为逆时针点旋转,R2为顺时针点旋转,R3为逆时针坐标系旋转,R4为顺时针坐标旋转
// float R1[2][2], R2[2][2], R3[2][2], R4[2][2];
// 计算 R1,R2,R3,R4, 用delta_theta;
get_rotate_matrix(...); // TODO( )
// 数值上,R3等于R2,R4等于R1
for (int i=0; i < iteration; i++) { // 迭代更新次数
//需要预测的点的个数集合:pts
// 移动点
Rotate(R1, pts);
Translate(pts, T);
// 旋转中心更新
rotate_center_update(pts);
// 移动坐标系
Rotate(R3, pts);
Translate(pts, T);
// 保存结果
save_projectory_prediction(iter, pts); // 保存结果时,需要反算点
}
// 运动学模型用来计算完Rt后,就可以Rt来一直迭代计算了
// /* x_r, y_r = rotate(R, point); 然后 t_x = dx - x_r; t_y = dy - y_r; 就是这样 */
// 坐标恢复
void coordinate_recover(...) {
for (int i=iter; i >= 0; i--) {
coordinate_system_recover(); // 顺时针Rt
rotate_center_recover(...); // 加上 减去的中心点坐标
}
rotate_center_recover(...); // 图像坐标系上 加上 减去的中心点坐标
}
void save_trajectory_point(int iter) {
coordinate_recover(iter, pts);
prediction_pts_result.push_back(pts); // 保存结果
}
void lateral_kinematics_model() {
// 用前/后轮模型都行,因为车辆上传感器可以获得,如下:
// 前后轮4个轮子的速度,前轮相对于车身的转角
// front_wheel_speed_left, front_wheel_speed_right
// rear_wheel_speed_left, rear_wheel_speed_right
// angle_for_front_wheel
// 我们取前轮两个轮子速度的平均值,作为单车模型的前轮速度,后轮同理。
// 以下用前轮举例
// 我们选取前轮作为模型的分析点,此时,有,如下模型公式
delta_theta = speed * sin(angle - M_PI / 2) / wheelbase * delta_time;
delta_x = speed * cos(angle + delta_theta) * delta_time;
delta_y = speed * sin(angle + delta_theta) * delta_time;
// 需要说明的是,angle换算成以水平x轴为起点旋转到当前方向的度数,
// 就是,相对于X轴的夹角。在模型里,是个相对的夹角,需要减去90度。
// delta_time 任取,取小点,比如 5ms,然后迭代算。
}
简单的方法:
伪代码,示意即可
void trajectory_update() {
// 每帧 结果清零 TODO ( )
// 每帧 模型的累加值清零
delta_theta = 0.f, delta_x = 0.f, delta_y = 0.f;
delta_time = _point_sparse_ratio_avm * delta_time_ori; // 系数超参数
int iter_num = int(delta_time / delta_time_ori);
// 要计算的车坐标点:4个轮胎,和车体边缘(车后视镜点)
std::vector<Point3D> points(m_vehicle_coords);
std::vector<Point3D> transformed_pts(m_vehicle_coords);
// 坐标中心换算,以车中心为原点
transform_coordinate_system(...);
// 这个 m_num_points_per_line 就是要预测的点的个数,预先给出
// 开始迭代
for (int i = 1; i <= num_points_per_line; i++) {
if (which_view == WhichView::kWide) {
// 如果是 去畸变 扳正后的正视图,这里是稀疏化点
iter_num = adjust_points_gap_ratio(i, delta_time);
}
// 这里就是迭代 N 次,一步到位,而不是 用大的delta_time(会error)
for (int j = 0; j < iter_num; j++) {
lateral_kinematics_model(delta_time_ori);
}
rotation_matrix(...); // 旋转矩阵
trans_vector(...); // 平移向量
// points清零,并用transformed_pts赋值 TODO ( )
for (int j = 0; j < 要计算的坐标点个数; j++) {
// 要计算的车坐标点:4个轮胎,和车体边缘(车后视镜点)
// R t
rotate(points[j], points[j]);
points[j].x += translation_vec[0];
points[j].y += translation_vec[1];
// 从车中心为原点 换算回图像坐标系
transform_coordinate_system_invert(...);
}
// 在这里保存点
}
}
极致的方法:
同心圆的代码就不贴了,实现起来还是很简单的。 就是要思路简单!实现也简单!
另外,关于刻度线 1米2米3米这种,用同心圆的方法计算即可。精确解
3.2 点的视图间转换
复杂的方法:
示意代码,伪代码
cv::Mat map_x(h, w, CV_32F);
cv::Mat map_y(h, w, CV_32F);
cv::Vec2f xy1, xy2;
int x_label = 0, y_label = 0;
for (int i = y_start; i < y_end; i++) {
float *data_x = map_x.ptr<float>(y_label);
float *data_y = map_y.ptr<float>(y_label);
for (int j = x_start; j < x_end; j++) {
矩阵乘法(H_bev_to_undistorted, j - x_start, i - y_start, xy1);
矩阵乘法(H_undistorted_to_thisview, xy1[0], xy1[1], xy2);
data_x[x_label] = xy2[0] - dw_bev_to_thisview;
data_y[x_label] = xy2[1] - dh_bev_to_thisview;
x_label++;
}
x_label = 0;
y_label++;
} // 得到bev to thisview的map:data_x, data_y
// 然后把 data_x, data_y 保存为bin
简单的方法:
H_bev_to_undistorted 直接和原始坐标coord0相乘,得到coord1,
然后 H_undistorted_to_thisview和coord1相乘,得到coord2,
coord2即为所求。
当在复用AVM标定流程的时候,直接调Homograpy即可,算都无需算。
而怎么计算 Homograpy单应性矩阵 这个,超出本文范围。如果掌握,请参看上面2.2给出的参考的链接
一些细节问题,还得实际中去慢慢踩坑,解决
3.3 点的稀疏化
在AVM的鸟瞰图BEV上,因为可以把点预测为均匀分布,而在 去畸变后的扳正后的图上,就难以做到均匀分布了。因为是拿BEV上的均匀分布的点 做Homograpy 透视变换 过去的,目前没对这个深究,根据肉眼观测,手动调节 iter_num 的大小,分段调节,比如前1个点,一次走多点,第2个点到第20个点一次走少点,远处的第100个点以上,再一次走多点,等等,超参数调节,把 透视变换后的不均匀 简单的处理了下,,这样使运算量较少的同时,使得 把点进行点与点之间连线,显得不太突兀,更平滑一点。当然,更好的方式,是拟合一个曲线,使画出来的点连线起来就是点密集分布平滑的线,因为可以指定Y坐标,逐渐增大,进行画点。
车辙轨迹预测,比起AVM整套来,还是太简单了,AVM流程一套弄下来,是永远的神.
以上写于2022年10月
void RutModel::rotate(Point3D &rotated_pt, Point3D origin_pt) {
float x = origin_pt.x;
float y = origin_pt.y;
rotated_pt.x = rotation_mat[0][0] * x + rotation_mat[0][1] * y;
rotated_pt.y = rotation_mat[1][0] * x + rotation_mat[1][1] * y;
}
void RutModel::lateral_kinematics_model(float delta_time) {
// rear wheel center mode
delta_theta +=
speed * tan(angle - m_pi / 2.f) / m_wheelbase_len * delta_time;
float r = m_wheelbase_len / tan(angle - m_pi / 2.f);
delta_x += (1.f - cos(delta_theta)) * r;
delta_y += sin(delta_theta) * r;
}
void RutModel::get_rotation_matrix() {
rotation_mat[0][0] = cos(delta_theta);
rotation_mat[0][1] = -1 * sin(delta_theta);
rotation_mat[1][0] = sin(delta_theta);
rotation_mat[1][1] = cos(delta_theta);
}
void RutModel::calc_trans_vector() {
Point3D object_point = m_rear_wheel_center;
object_point.x -= m_vehicle_center.x;
object_point.y -= m_vehicle_center.y;
object_point.y *= -1.f;
Point3D predicted_pt;
predicted_pt.x = object_point.x + delta_x;
predicted_pt.y = object_point.y + delta_y;
Point3D rotated_pt;
rotate(rotated_pt, object_point);
translation_vec[0] = predicted_pt.x - rotated_pt.x;
translation_vec[1] = predicted_pt.y - rotated_pt.y;
}
// 部分实现
...
get_rotation_matrix();
calc_trans_vector();
points.clear();
points.assign(transformed_pts.begin(), transformed_pts.end());
for (int j = 0; j < num_coords; j++) {
if (filter_flag[j]) {
points[j].x = m_invalid_value;
points[j].y = m_invalid_value;
continue;
}
rotate(points[j], points[j]);
points[j].x += translation_vec[0];
points[j].y += translation_vec[1];
points[j].y *= -1;
points[j].x += m_vehicle_center.x;
points[j].y += m_vehicle_center.y;
boundary_x[j] = points[j].x;
boundary_y[j] = points[j].y;
}
trajectory_points.push_back(points);
}