纯定位模式 (一)

纯定位时载入的地图在submap中的trajectory为0,启动定位后,系统默认从起点开始建立子图进行匹配,该子图在submap中的trajectory为 1。

一条trajectory可以理解为一次建图过程。如果只让机器人跑一圈,那么trajectory数就是1. 但是,建好图后又需要在图中走,这时候可以增加一条trajectory,把pure_localization设为true,那么机器人再重新跑的过程中就会跟已经建好的图进行匹配,估计机器人在地图中的路径。所以,一次运行就代表了一条trajectory。

cartographer定位时,它的地图是变化的

纯定位和建图不同,local和global SLAM之间的延迟要尽量低. 其次global SLAM通常会在作为地图的frozen trajectory和当前trajectory之间,找到大量的相互约束(inter constraints)。

纯定位模式默认只保存 3 个子图,只有这些子图内的节点参与优化。所以地图大小并不会影响后端优化的速度。而建图结束的时候全图优化应该是所有节点参与。

我们需要显著减小POSE_GRAPH.optimize_every_n_nodes以接收频繁的结果, global SLAM通常会太慢,所以再显著减小global_sampling_ratioconstraint_builder.sampling_ratio以补偿大量的约束。

submaps.resolution应当和运行的.pbstream中的子图resolution一致

建图时,cartographer使用 fast correlative scan matcher找到回环,search window默认为7m×7m×30°. 但纯定位模式中,search window变得非常大,所以纯定位消耗CPU更大

constraint_builder_2d.cc中的ComputeConstraint()函数,在纯定位模式下,会频繁地调用MaybeAddGlobalConstraint(),这个函数会搜索整个子图,造成的开销特别大。建议将global_sampling_ratio调小,这样就可以从一定程度上减少开销。

数据集的纯定位测试

程序运行时只在局部建图(保留固定数量的子图),随着机器人的运动和观测, 新的子图会添加进来,旧的子图会被删除,就是只保留局部的地图数据。我的理解是保留局部数据越多,与原有地图进行回环检测就增多,这样计算量增加的同时也增加了误匹配的概率,使得整个系统不稳定 ,所以这里只进行定位,不用于建图

需要用到数据集b2-2016-04-05-14-44-52.bagb2-2016-04-27-12-31-41.bag,前者用于建图,后者用于定位

先运行建图:

1
roslaunch cartographer_ros offline_backpack_2d.launch bag_filenames:=${HOME}/path/b2-2016-04-05-14-44-52.bag

等一会,完成后会生成pbstream文件,如果运行建图,实际地图是这样的
建图.png

这里是先用第一个数据集建图,建完图后用第二个数据集进行定位测试。

1
roslaunch cartographer_ros demo_backpack_2d_localization.launch load_state_filename:=${BAG}/b2-2016-04-05-14-44-52.bag.pbstream bag_filename:=${BAG}/b2-2016-04-27-12-31-41.bag

后面两个参数一个是上一步建好的地图,pbstream格式,另一个是数据集。你会看到机器人很快实现了定位

在机器人的旧版本纯定位

用于定位时的home_no_odom_localization.launch

1
2
3
4
5
6
7
8
9
10
11
12
<launch>
<node name="cartographer_node" pkg="cartographer_ros"
type="cartographer_node" args="
-configuration_directory $(find cartographer_ros)/configuration_files
-configuration_basename home_no_odom_localization.lua
-load_state_filename /home/xiaoqiang/Documents/ros/maps/home.pbstream"
output="screen">
</node>

<node name="cartographer_occupancy_grid_node" pkg="cartographer_ros"
type="cartographer_occupancy_grid_node" args="-resolution 0.05" />
</launch>

使用cartographer建图和定位时,launch文件的不同之处在于cartographer_node节点的参数-configuration_basename-load_state_filename,前者就是lua文件的不同,后者只在定位时才有,毕竟没有地图就不能定位

建图使用no_odom_and_imu.lua,定位使用home_no_odom_localization.lua:

1
2
3
4
5
6
7
include "home.lua"

TRAJECTORY_BUILDER.pure_localization = true
--- 每10个有效帧一个子图,子图构建完成要闭环检测一次,这个数越小,闭环检测越频繁,当然CPU增大
POSE_GRAPH.optimize_every_n_nodes = 5
POSE_GRAPH.constraint_builder.sampling_ratio = 0.1
POSE_GRAPH.global_sampling_ratio = 0.002

启动roslaunch cartographer_ros home_localization.launch后应该就可以了,必须有map坐标系,否则就是没成功。
定位结果.png

optimize_every_n_nodes会影响定位时的CPU占用,设置为5时对应50%,但只有定位时的一秒左右,所以无所谓。

添加轨迹实际调用的源码是node.cc中的Node::HandleStartTrajectory——AddTrajectory

新版本的配置

pure_localization_trimmer.png
这是因为我的 cartographer版本旧,不支持pure_localization_trimmer,最好还是更新到最新。在18年7月增加了pure_localization_trimmer,删掉了pure_localization,如果是新一点的版本,做对应修改。commit历史在这里

相应的源码在 MapBuilder::AddTrajectoryBuilder——MaybeAddPureLocalizationTrimmer

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
if (trajectory_options.pure_localization()) 
{
LOG(WARNING) << "'TrajectoryBuilderOptions::pure_localization' field is deprecated"
"Use 'TrajectoryBuilderOptions::pure_localization_trimmer' instead.";
pose_graph->AddTrimmer(absl::make_unique<PureLocalizationTrimmer>(
trajectory_id, 3 /* max_submaps_to_keep */) );
return;
}
if (trajectory_options.has_pure_localization_trimmer())
{
// 添加一个PureLocalizationTrimmer类型修饰器
pose_graph->AddTrimmer(absl::make_unique<PureLocalizationTrimmer>(
trajectory_id,
trajectory_options.pure_localization_trimmer().max_submaps_to_keep() ) );
}

新版本的配置:

1
2
3
4
TRAJECTORY_BUILDER.pure_localization_trimmer = {
--- 最大保存子图数,存定位模式通过子图进行定位,但只需要当前和上一个子图即可
max_submaps_to_keep = 3,
}

如果使用了机器人的odom,但无法定位成功或定位误差太大,可能的原因是odom不够精确,尝试关闭odom,并由cartographer代为提供odom。既然定位误差大,可能在建图时的原始submap之间误差较大,图优化后map的边缘较为模糊

重定位注释掉Node::GrawAndPublish函数的后三行,否则会继续建图。但是建图的时候还得加回来。
由于Cartographer定位是当前帧和子图进行匹配,所以对于机器人周围被新物体遮挡的情况完全没有影响。

唯一影响就是没法重定位了,因为地图上没有这些新物体。但是就算长时间被遮挡,存在累积定位误差,只要机器人离开了遮挡区域,回到了之前的环境,就可以重定位回来。

获得机器人在map坐标系中的坐标

cartographer是使用tf转换来发布位姿的,如果想发布类似amcl_pose那样的话题,只能修改源码或自己写个程序。

修改源码的方法可以参考这里

我还是自己写了个程序,其实就是tf报警:Lookup would require extrapolation into the past中的程序改了坐标系和话题名称,以获得base_footprintmap中的坐标系

问题

某次运行得到这样结果:

追查到node.cc里的Node::FinishTrajectoryUnderLock函数,检查了一下是否可以关掉,指定id是否存在,是否已经被Finished了等情况后,如果一切正常,则停止订阅Topic、清除id及其他与该trajectory相关的量。最后调用map_builder_bridge_中的FinishTrajectory函数。发现code对应cartographer_ros_msgs::StatusCode::INVALID_ARGUMENT

参考:
官方参考
更优雅的设定初始位姿


Realsense D435i关闭IR结构光

由于要做Realsense D435i的双目结构光相机标定,其中用到了ROS来录制数据包,但是结构光会影响标定,所以得先关闭IR结构光发射器。

一次性关闭IR光

其实就是动态参数调整

1
2
roslaunch realsense2_camera rs_camera.launch 
rosrun rqt_reconfigure rqt_reconfigure

rqt_reconfigure
emitter_enabled有三项: Laser(1), Laser Auto(2), Off(0)。 选Off即可

通过肉眼看相机发射器已经不再发射IR光了,打开rqt,找到/camera/infra1/image_rect_raw/camera/infra2/image_rect_raw,发现确实没有白斑了。

永久关闭

其实就是通过realsense-viewer,这个是驱动层面的操作,更加底层
配置栏.png

参考:Realsense D435i 关闭IR光


ROS的topic_tools包

throttle

转播 (relay)话题, 限制最大发布频率或者带宽。使用格式: rosrun topic_tools throttle messages <intopic> <msgs_per_sec> [outtopic]

messages是当前使用的是throttle的messages模式,限制发布频率. intopic是指你想要改变频率的那个topic, msgs_per_sec是指你想要它发布的频率,而outtopic是指改变发布频率后的topic的名称,可以省略,如果省略则自动在原来topic的名字上后缀throttle

另外还有bytes模式, 用以限制带宽: rosrun topic_tools throttle bytes <intopic> <bytes_per_sec> <window> [outtopic]

例如,让雷达的带宽占用降至1KBps,则命令为:

1
rosrun topic_tools throttle bytes base_scan 1024 1.0

改变topic发布频率并不是说原来的topic就没了,或者直接在原来的topic上更改, 而是把其更改后的topic发布出来,原来的topic仍然存在

有三个参数,需要注意的是~lazy,如果它等于True的话,只有当你订阅到throttle发出来的消息,它才会工作,这显然是ros::publish函数的lazy模式了

最大的问题是 频率控制的精度低,我要求4Hz,实际却在3.6Hz左右,所以实际要设置的大一些

transform

订阅一个topic或topic field,并在应用给定的Python表达式后将传入的数据发布到另一个topic。它可以处理任何消息类型,主要用于简单的消息转换,例如计算向量或四元数的范数, 甚至将四元数转换为欧拉角。

transform <input> <output_topic> <output_type> [<expression>] [--import <module> [<module> ...]]

  • input:要订阅的传入topic(或topic field)
  • output_topic:要发布的输出topic
  • expression:转换输入消息的Python表达式,在变量m中给出。默认表达式是m,它将输入(可以是topic field)转发到output_topic
  • import :要在表达式中导入和使用的Python模块的列表。默认导入numpy模块
1
2
3
4
# 计算IMU给出的方向四元数的范数
rosrun topic_tools transform /imu/orientation /norm std_msgs/Float64 'numpy.linalg.norm([m.x, m.y, m.z, m.w])'
# 将方向四元数转换为Euler角度
rosrun topic_tools transform /imu/orientation /euler geometry_msgs/Vector3 'tf.transformations.euler_from_quaternion([m.x, m.y, m.z, m.w])' --import tf


  • mux: multiplex between multiple topics

  • relay: republish data on one topic to another.

  • relay_field: allow to republish data in a different message type

  • drop: relay a topic, dropping X out of every Y message.

mux

topic_tools-mux
可以用于多地图实时切换,参考MUX的使用

参考:throttle Wiki


标定D435i(二) 相机

Kalibr可以解决以下的标定问题:

  • 多相机标定: 一个相机系统的内外参标定,这几个相机没有全局性重叠的视角
  • 视觉惯性标定(camera-IMU): IMU关于相机系统的时空间标定
  • Rolling Shutter Camera calibration: full intrinsic calibration (projection, distortion and shutter parameters) of rolling shutter cameras

可以下载Kalibr源码编译生成可执行文件,也可以下载其CDE精简版包。这中间有个坑就是CDE精简包是没有办法标定彩色图片的,而D435输出的是彩色图。所以还是按编译源码的方式

使用Kalibr标定相机的内参和多个相机相对位置关系即外参

安装Kalibr

  1. 安装依赖项

    1
    sudo apt-get install python-pyx python-setuptools python-rosinstall ipython libeigen3-dev libboost-all-dev doxygen ros-kinetic-vision-opencv ros-kinetic-image-transport-plugins ros-kinetic-cmake-modules python-software-properties software-properties-common libpoco-dev python-matplotlib python-scipy python-git python-pip ipython libtbb-dev libblas-dev liblapack-dev python-catkin-tools libv4l-dev
  2. 源码放入工作空间进行编译,会花很长时间,所以编译命令要这样:

    1
    2
    # 视情况取j8
    catkin_make -DCMAKE_BUILD_TYPE=Release -j4

编译kalibr可能会出现 fatal error: numpy/arrayobject.h: No such file or directory ,解决方法: sudo apt-get install --reinstall python-numpy

可能出现catkin_make时,下载suitesparse过久甚至失败的问题。解决方法: 修改~/catkin_ws/src/kalibr/suitesparse中的CMakeLists.txt为新CMakeLists.txt, 然后重新catkin_make

标定板

下载标定板

我下载的是Aprilgrid 6x6 0.5x0.5m(unscaled), 打印在A3纸上

原版的参数是:6X6 tags,6乘6个格子,一个大格子size=5.5cm,一个小格子spacing=1.65cm

A3纸上的缩放:6X6 tags,一个大格子size=3.5cm,一个小格子spacing=1cm。记得打印出来用尺子量一下,以免出现差错。
标定板参数示意图

下载官网提供的yaml格式文件,需要按照设定的尺寸进行修改

1
2
3
4
5
target_type: 'aprilgrid' #gridtype
tagCols: 6 #number of apriltags
tagRows: 6 #number of apriltags
tagSize: 0.035 #size of apriltag, edge to edge [m]
tagSpacing: 0.3 #ratio of space between tags to tagSize

找一个适合的能拍到棋盘格的距离,启动相机: roslaunch realsense2_camera rs_camera.launch

d453i是有红外发射器的,可以发射很多红外小斑点,如果打开你会在rviz看到很多光斑,可能不利于标定,所以标定时关闭这个发射器的。

降低图像话题的频率,录制图像数据包

kalibr在处理标定数据的时候要求图像的频率不可过高,一般为4hz(后面计算过程报错,改为20)。使用如下指令来限制图像频率:

1
2
3
rosrun topic_tools throttle messages /camera/color/image_raw 4 /color
rosrun topic_tools throttle messages /camera/infra2/image_rect_raw 4 /infra_left
rosrun topic_tools throttle messages /camera/infra1/image_rect_raw 4 /infra_right

用topic_tools throttle限制频率后,一定要查看限制后的topic输出频率:rostopic hz /topic,你会发现实际的频率与设定的频率并不一致,比如:rosrun topic_tools throttle messages /topic_1 25 /topic_2,如果topic_1是40hz,/topic_2可能不是25hz,而是20hz,而且每次的实际频率可能不同,具体原因不明。

注意这里是采用了新的话题/color去发布,所以下面录制要写/color话题

将标定目标AprilGrid置于相机前方合理距离范围内,准备录制两分钟:

1
rosbag record -O multicameras_calibration /infra_left /infra_right /color

开始缓慢移动标定板,让所有摄像头看到标定板不要太远,不然无法检测到标定目标的特征,这个过程中看不到图案的识别结果,这是kalibr的原因。在标定算法中需要检测是否有足够数量图片检测到标定特征,否则直接无法标定。移动标定物时候不要过快导致运动模糊,我们只需要获取不同位置和角度的图像,确保图像清晰和特征完整即可。另外要尽可能多角度和多位置(上下左右等)甚至到摄像头捕捉图像的边缘,这样移动目标1min左右即可。

移动的方式可以参考TUM CALIBRATION SEQUENCES的标定方式,点play即可播放

kalibr算法计算各个摄像头的内参和外参

april_6x6_A3.yaml: 标定物的参数,具体是标定目标的尺寸之类,因为我是缩小打印在A3上,所以要对参数进行修改;pinhole-equi – 选择的相机模型,kalibr提供了很多相机模型,可以自己选择; —bag-from-to 可以选择时间段,毕竟录制的时候不能保证整体都录制的很好。标定会花很长时间,最后会输出一个pdf和txt文件,有内外参数据。

只标定主相机:

1
kalibr_calibrate_cameras --target ../yaml/april_6x6_A4.yaml --bag ./bag/0_multicameras_calibration.bag --model pinhole-equi  --topic  /color  --show-extraction --approx-sync 0.04

最后还是标定的多相机:

1
kalibr_calibrate_cameras --target april_6x6_A3.yaml --bag multicameras_calibration.bag --models pinhole-equi pinhole-equi pinhole-equi --topics /infra_left /infra_right /color  --show-extraction --approx-sync 0.04 --bag-from-to 10 100

  • pinhole-radtan指的是针孔相机模型和畸变模型,每个相机都要指定。还有Pinhole + FOV等等
  • --bag-from-to 10 100指的是录制的第26秒到100秒这段时间
  • --show-extraction, 在标定过程中可视化角点检测情况是否良好
  • --approx-sync 0.04

结果报错: ImportError: No module named igraph
解决方法: sudo apt-get install python2.7-igraph

需要在有界面的情况下标定,因为会弹出几个窗口,所以不能通过SSH进行,类似这样:

可以使用calibration validator进行标定的验证,原理是对重投影误差进行量化分析,同样需要有界面:
kalibr_camera_validator --cam camchain-multicameras_calibration.yaml --target april_6x6_A3.yaml
验证结果.png


标定D435i(一) IMU

realsense d435i包含两个红外相机、红外发射器、RGB相机和IMU四个模块,显然四个传感器的空间位置是不同的,我们在处理图像和IMU数据时需要将这些数据都放在统一的坐标系上去。比如我们用d435i运行vins,处理的图像和IMU数据都需要放在同一个坐标系下,因此需要标定IMU相对RGB相机的空间位置(包括旋转和位移)。

另外,相机固有参数比如焦距、畸变参数等以及IMU的零偏和scale系数等都需要提前知道。前者称为外参,后者称为内参,在运行程序前我们需要标定它们,不论程序是否有自标定功能,毕竟好的初始标定值对于自标定来说也是有利的。

标定顺序:IMU标定 —> 相机标定 —> IMU+相机联合标定. 这么设定顺序是因为最后一步的IMU和相机的联合标定需要 IMU和相机的内参

Allan方差

在IMU采集数据时,会产生两种误差:确定性误差和随机性误差,为获得精确的数据,需要对上述两种误差进行标定,加速度计和陀螺仪随机误差的标定通常使用Allan方差法,它是20世纪60年代由美国国家标准局的David Allan提出的基于时域的分析方法。

A ROS package tool to analyze the IMU performance. C++ version of Allan Variance Tool. The figures are drawn by Matlab, in scripts.

Actually, just analyze the Allan Variance for the IMU data. Collect the data while the IMU is Stationary, with a two hours duration.
code_utils标定IMU的噪音密度和随机游走系数。

安装库和依赖项

安装依赖项,不装之后的编译会报错: sudo apt-get -y install libdw-dev,结果可能提示

1
2
3
4
The following packages have unmet dependencies:
libdw-dev : Depends: libelf-dev but it is not going to be installed
Depends: libdw1 (= 0.165-3ubuntu1) but it is not going to be installed
E: Unable to correct problems, you have held broken packages.

这是因为一个依赖项已经安装了不同版本:Depends: libelf1 (= 0.165-3ubuntu1) but 0.165-3ubuntu1.2 is to be installed。 解决方法: sudo aptitude install libdw-dev,对给出的方案,选择第二个,降级 libelf1[0.165-3ubuntu1.1 (now) -> 0.158-0ubuntu]

编译code_utils和imu_utils

全局安装ceres库,因为code_imu依赖ceres。不要同时把imu_utils和code_utils一起放到src下进行编译。因为imu_utils 依赖 code_utils,原作者的CMakeLists写的不好,所以先编译code_utils再编译后者。

code_utils下面找到sumpixel_test.cpp,修改#include "backward.hpp"#include"code_utils/backward.hpp",再依次编译两个包

发布D435i的IMU数据

可以直接在rs_camera.launch基础上针对IMU校准做修改。目的是将acc、gyro数据对齐使用同一个topic发布。

1
2
3
4
<!-- 更改前原版本arg name="unite_imu_method"          default=""/-->  
<arg name="unite_imu_method" default="linear_interpolation"/>
<!--或着将参数改为copy-->
<arg name="unite_imu_method" default="copy"/>

启动: roslaunch realsense2_camera rs_imu_calibration.launch,然后录制imu数据包rosbag record -O imu_calibration /camera/imu,让IMU静止不动两个小时,录制IMU的bag.

标定

根据imu_utils文件夹里面的A3.launch改写D435i标定启动文件:d435i_imu_calib.launch注意,记得修改max_time_min对应的参数,默认是120,也就是两个小时,如果ros包里的imu数据长度没有两个小时,等bag播放完了,还是停留在wait for imu data这里,不会生成标定文件。我录了1小时59分多一点,所以还得改成119

d435i_imu_calibration.launch如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
<launch>
<node pkg="imu_utils" type="imu_an" name="imu_an" output="screen">
<!--TOPIC名称和上面一致-->
<param name="imu_topic" type="string" value= "/camera/imu"/>
<!--imu_name 无所谓-->
<param name="imu_name" type="string" value= "D435i"/>
<!--标定结果存放路径-->
<param name="data_save_path" type="string" value= "$(find imu_utils)/data/"/>
<!--数据录制时间-min-->
<param name="max_time_min" type="int" value= "120"/>
<!--采样频率,即是IMU频率,采样频率可以使用rostopic hz /camera/imu查看,设置为200,也就是rosbag play播放频率-->
<param name="max_cluster" type="int" value= "200"/>
</node>
</launch>

先启动标定程序: roslaunch imu_utils d435i_imu_calib.launch,再播放bag: rosbag play -r 200 imu_calibration.bag

标定结果

标定结果是D435i_imu_param.yaml:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
%YAML:1.0
---
type: IMU
name: D435i
Gyr:
unit: " rad/s"
avg-axis:
gyr_n: 3.6673681012835031e-03
gyr_w: 7.0017785520472972e-05
x-axis:
gyr_n: 3.6001489799186333e-03
gyr_w: 6.2846247607788020e-05
y-axis:
gyr_n: 4.7157261366663813e-03
gyr_w: 7.5207268006344615e-05
z-axis:
gyr_n: 2.6862291872654953e-03
gyr_w: 7.1999840947286307e-05
Acc:
unit: " m/s^2"
avg-axis:
acc_n: 2.7436489578044256e-02
acc_w: 1.0915021608117670e-03
x-axis:
acc_n: 1.8271976632141730e-02
acc_w: 5.5394830052109354e-04
y-axis:
acc_n: 2.8924134998445018e-02
acc_w: 1.5674764920646303e-03
z-axis:
acc_n: 3.5113357103546017e-02
acc_w: 1.1530816898495772e-03

我们一会只用到Gyr中的avg-axisgyr_ngyr_w, Acc中的avg-axisacc_nacc_w

终端输出结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
gyr x  numData 781205
gyr x start_t 1.59343e+09
gyr x end_t 1.59344e+09
gyr x dt
-------------7140.59 s
-------------119.01 min
-------------1.9835 h
gyr x freq 109.403
gyr x period 0.00914049
gyr y numData 781205
gyr y start_t 1.59343e+09
gyr y end_t 1.59344e+09
gyr y dt
-------------7140.59 s
-------------119.01 min
-------------1.9835 h
gyr y freq 109.403
gyr y period 0.00914049
gyr z numData 781205
gyr z start_t 1.59343e+09
gyr z end_t 1.59344e+09
gyr z dt
-------------7140.59 s
-------------119.01 min
-------------1.9835 h
gyr z freq 109.403
gyr z period 0.00914049
Gyro X
C -6.83161 94.2973 -19.0588 2.983 -0.0404918
Bias Instability 2.37767e-05 rad/s
Bias Instability 6.28462e-05 rad/s, at 63.1334 s
White Noise 12.9453 rad/s
White Noise 0.00360015 rad/s
bias -0.363298 degree/s
-------------------
Gyro y
C -8.74367 117.584 -15.9277 2.47408 -0.0373467
Bias Instability 6.41864e-05 rad/s
Bias Instability 7.52073e-05 rad/s, at 104.256 s
White Noise 16.8998 rad/s
White Noise 0.00471573 rad/s
bias -0.544767 degree/s
-------------------
Gyro z
C -4.51808 68.1919 -9.33284 1.95333 -0.0262641
Bias Instability 8.50869e-05 rad/s
Bias Instability 7.19998e-05 rad/s, at 63.1334 s
White Noise 9.43212 rad/s
White Noise 0.00268623 rad/s
bias -0.0762471 degree/s
-------------------
==============================================
==============================================
acc x numData 781205
acc x start_t 1.59343e+09
acc x end_t 1.59344e+09
acc x dt
-------------7140.59 s
-------------119.01 min
-------------1.9835 h
acc x freq 109.403
acc x period 0.00914049
acc y numData 781205
acc y start_t 1.59343e+09
acc y end_t 1.59344e+09
acc y dt
-------------7140.59 s
-------------119.01 min
-------------1.9835 h
acc y freq 109.403
acc y period 0.00914049
acc z numData 781205
acc z start_t 1.59343e+09
acc z end_t 1.59344e+09
acc z dt
-------------7140.59 s
-------------119.01 min
-------------1.9835 h
acc z freq 109.403
acc z period 0.00914049
acc X
C 3.36177e-05 0.00175435 -0.000159698 7.23303e-05 -7.16006e-07
Bias Instability 0.000553948 m/s^2
White Noise 0.018272 m/s^2
-------------------
acc y
C 9.36955e-05 0.00234733 0.00012197 0.000243676 -2.66252e-06
Bias Instability 0.00156748 m/s^2
White Noise 0.0289241 m/s^2
-------------------
acc z
C 5.07832e-05 0.00331104 -0.000381222 0.000199602 -2.43776e-06
Bias Instability 0.00115308 m/s^2
White Noise 0.0351134 m/s^2

经过这些标定会生成一个yaml文件和很多txt文件,主要是yaml文件,给出了加速度计和陀螺仪三轴的noise_densityrandom_walk,同时计算出了平均值,后面IMU+摄像头联合标定的时候需要这些均值。

标定外参的准备

将Acc和Gyr的第一组平均数据拷贝到kalibr对应的imu.yaml文件中

1
2
3
4
5
6
7
rostopic: /camera/imu
update_rate: 200.0 #Hz

accelerometer_noise_density: 2.89e-01 #continous
accelerometer_random_walk: 4.55e-04
gyroscope_noise_density: 3.02e-03 #continous
gyroscope_random_walk: 2.29e-05

分别是加速度计和陀螺仪的高斯白噪声和随机游走的平均值,是IMU噪声模型中的两种噪声。

 查看默认imu与相机参数

D435i相关的imu_info话题如下:
imu_info列表

/camera/gyro/imu_info和/camera/accel/imu_info没有结果

realsense_camera/IMUInfo with the header.frame_id set to either imu_accel or imu_gyro to distinguish between accel and gyro info. 消息成员如下:

1
2
3
4
string frame_id
float64[12] data
float64[3] noise_variances
float64[3] bias_variances

rostopic echo -n1 /camera/color/camera_info得到结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
  frame_id: "camera_color_optical_frame"
height: 480
width: 640
distortion_model: "plumb_bob"
D: [0.0, 0.0, 0.0, 0.0, 0.0]
K: [611.3538208007812, 0.0, 327.437744140625, 0.0, 610.015869140625, 239.99667358398438, 0.0, 0.0, 1.0]
R: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
P: [611.3538208007812, 0.0, 327.437744140625, 0.0, 0.0, 610.015869140625, 239.99667358398438, 0.0, 0.0, 0.0, 1.0, 0.0]
binning_x: 0
binning_y: 0
roi:
x_offset: 0
y_offset: 0
height: 0
width: 0
do_rectify: False

对于/camera/depth/imu_info,结果是:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
frame_id: "camera_depth_optical_frame"
height: 480
width: 640
distortion_model: "plumb_bob"
D: [0.0, 0.0, 0.0, 0.0, 0.0]
K: [380.23321533203125, 0.0, 316.4999084472656, 0.0, 380.23321533203125, 237.40985107421875, 0.0, 0.0, 1.0]
R: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
P: [380.23321533203125, 0.0, 316.4999084472656, 0.0, 0.0, 380.23321533203125, 237.40985107421875, 0.0, 0.0, 0.0, 1.0, 0.0]
binning_x: 0
binning_y: 0
roi:
x_offset: 0
y_offset: 0
height: 0
width: 0
do_rectify: False

话题/camera/infra2/camera_info/camera/infra2/camera_info的结果完全一样:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
frame_id: "camera_infra1_optical_frame"
height: 480
width: 640
distortion_model: "plumb_bob"
D: [0.0, 0.0, 0.0, 0.0, 0.0]
K: [380.23321533203125, 0.0, 316.4999084472656, 0.0, 380.23321533203125, 237.40985107421875, 0.0, 0.0, 1.0]
R: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
P: [380.23321533203125, 0.0, 316.4999084472656, 0.0, 0.0, 380.23321533203125, 237.40985107421875, 0.0, 0.0, 0.0, 1.0, 0.0]
binning_x: 0
binning_y: 0
roi:
x_offset: 0
y_offset: 0
height: 0
width: 0
do_rectify: False

可以看出,实际上出厂没有标定IMU

绘制Allan曲线

imu_utils/scripts中是Matlab写的.m文件,按照draw_allan.m创建文件draw_allan_acc.mdraw_allan_gyr.m

由于Ubuntu下安装Matlab比较麻烦,因此将数据和.m文件都拷贝到Windows系统下绘制,注意文件路径和m文件中要对应,尤其matlab中的当前路径指的是左侧栏的路径,而不是m文件所在的路径

allan_gyr.jpg
allan_acc.jpg

参考:
RealSense D435i Calibration
官方Github
VIO标定IMU随机误差——Allan方差法


标定D435i(三) 外参数
  • 标定板pattern尽量选择apriltag,尽量不要用纸质的
  • 选择合适的相机模型
  • 标定手法看看TUM的

realsense D435i安装配置


realsense D435i终于到手了,打开发现其实很小巧,先做一些配置看看。默认要使用USB3,如果RealSense使用USB2, Output Resolution自动降到480*270 30fbs,而非产品所宣称的1280 x 720 active stereo depth resolution和90fps,且只有Stereo Moudle在工作,无Image Sensor的RGB Moudle菜单项,无法进行3D建模。

D435i采用了即用型USB供电形式,不仅提供深度传感器模组,还配备了一个IMU单元(惯性测量单元,采用的博世BMI055)。凭借内置的IMU单元,结合视觉数据可实现6DoF追踪功能。其中,IMU将各种线性加速度计和陀螺仪数据结合,可检测X,Y,Z三轴的旋转和平移,以及俯仰、横摇等动作。D435i的2000万像素RGB摄像头和3D传感器可以30帧/秒的速度提供分辨率高达1280 × 720,或者以90帧/秒的速度提供848 × 480的较低分辨率。该摄像头具有全局快门,可以处理快速移动物体,室内室外皆可操作。深度距离在0.1 m~10 m之间,视场角度为85 × 58度

可以获得RGB图、左右红外摄像图、深度图、IMU数据,并且将深度图数据和RGB图进行对齐。左右红外相机进行测量深度,中间红外点阵投射器相当于补光灯,不打开也能测深度,只是效果不好;最右边的rgb相机用于采集彩色图片,最终可以将彩色视频流与深度流进行对齐.

两个驱动讲究搭配,否则安装编译都不报错,运行时会报错 realsense hwmon command 0x7d failed 我使用的是realsense-ros-2.2.23librealsense-2.43.0

安装librealsense

先装realsense的驱动,步骤参考网上的,其实不必完全相同,我的步骤如下:

  1. 下载驱动,然后解压到根目录
  2. 执行如下命令:
    1
    2
    3
    4
    5
    cd ~/librealsense
    sudo apt-get install libssl-dev libusb-1.0-0-dev pkg-config libgtk-3-dev libglfw3-dev
    # 许可脚本
    sudo cp config/99-realsense-libusb.rules /etc/udev/rules.d/
    sudo udevadm control --reload-rules && udevadm trigger
  3. 开始编译
    1
    2
    3
    4
    5
    cd ~/librealsense
    mkdir build && cd build
    cmake ../ -DBUILD_EXAMPLES=true

    sudo make uninstall && make clean && make -j7 && sudo make install

运行realsense-viewer验证,看到如下画面
RGB Depth Infrared
配置可以修改为常见的黑白深度图以及分辨率等等
配置栏

ROS驱动

这里下载,然后放到某工作空间,编译即可:

1
2
catkin_make -DCATKIN_ENABLE_TESTING=False -DCMAKE_BUILD_TYPE=Release
catkin_make install

现在运行roslaunch realsense2_camera rs_camera.launch,打开rqt就能看到realsense的三种图像
实际运行了/camera/realsense2_cameracamera/realsense2_camera_manager两个节点,涉及的话题有很多,以后慢慢分析。

tf关系如图:
frames.png

在rqt中打开深度图时出现了报错,又是图像编码问题
launch文件中报错.png
rqt终端报错.png
在rqt里查看/camera/depth/image_rect_raw/theora,同样没有深度图,但rqt不报警

在rqt查看RGB和红外的图像,只要选择theora,rqt终端都会报警:

1
[ WARN][/rqt_gui_cpp_node_25764] [TheoraSubscriber::internalCallback] line_170  [theora] Packet was not a Theora header


使用对齐的深度话题信息发布RGBD点云

1
roslaunch realsense2_camera rs_rgbd.launch

这一步还没有成功,在rviz里没看到结果

内参

realsense启动时可以发现信息 [color stream is enabled - width: 640, height: 480, fps: 15, Format: RGB8]

realsense相机出厂的时候一般都标定好了,直接读取他们的内参即可。终端输入: rs-sensor-control,根据提示选择,出现所有视频流的列表时,根据上面的信息选择,最终显示内参:

1
2
3
4
Principal Point         : 327.438, 239.997
Focal Length : 611.354, 610.016
Distortion Model : Inverse Brown Conrady
Distortion Coefficients : [0,0,0,0,0]

疑难问题 undefined symbol: _ZN2cv3MatC1Ev

运行rs_camera.launch时报错

1
symbol lookup error: /home/jetson/catkin_ws/devel/lib//librealsense2_camera.so: undefined symbol: _ZN2cv3MatC1Ev   [camera/realsense2_camera_manager-2] process has died 

我查了好长时间,期间还改用其他的librealsense和firmware的版本,比如到firmware版本历史,下载对应版本。 参考firmware-update-tool的Usage部分

运行时可以看到相关的环境版本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
RealSense ROS v2.3.2
[ INFO] [1723892537.013079361]: Built with LibRealSense v2.50.0
[ INFO] [1723892537.013109698]: Running with LibRealSense v2.50.0
[ INFO] [1723892537.311113159]: Device with serial number 014122072296 was found.
[ INFO] [1723892537.311624883]: Device with physical ID 2-1.4-6 was found.
[ INFO] [1723892537.311852281]: Device with name Intel RealSense D435I was found.
[ INFO] [1723892537.313360062]: Device with port number 2-1.4 was found.
[ INFO] [1723892537.313859786]: Device USB type: 3.2

[ INFO] [1723892537.373681930]: JSON file is not provided
[ INFO] [1723892537.373700779]: ROS Node Namespace: camera
[ INFO] [1723892537.373737260]: Device Name: Intel RealSense D435I
[ INFO] [1723892537.373757100]: Device Serial No: 014122072296
[ INFO] [1723892537.373777005]: Device physical port: 2-1.4-6
[ INFO] [1723892537.373792685]: Device FW version: 05.13.00.50

后来又参考如何修改CMakeLists

最后发现其实存在两个librealsense2_camera.so

1
2
/home/jetson/catkin_ws/devel/lib/librealsense2_camera.so
/opt/ros/noetic/lib/librealsense2_camera.so

只需要运行ROS自带的就可以了: roslaunch /opt/ros/noetic/share/realsense2_camera/launch/rs_camera.launch

参考:
realsense D435 驱动安装并嵌入ROS中使用


(九) 总结

gmapping源码中不完善的部分:

  1. drawFromMotion函数是一个十分粗糙的运动模型,只是简单的矢量加减运算。相比于《概率机器人》中提到的速度模型和里程计模型,有很多方面都没有考虑,精度上可能有折扣。
  2. m_angularOdometryReliabilitym_linearOdometryReliability用于控制激光评分的增益odo_gain, 实际都是0,所以没有用到,而且无法用参数赋值, 要使用就得修改源码.
  3. 紧接上面的代码,localScore=odo_gain*score(map, localPose, readings);,localScore在之后的likelihoodAndScore函数里实际赋值为0,根本没有用。所以这里有好几行可以删除
  4. 在论文中,粒子的权重不是用最优位姿的似然值来表示的,而是用所有的似然值的和来表示的
  5. 重采样的粒子索引选取用的是轮盘赌算法,有些论文提出了一些更好的算法

OpenMP加速gmapping

在CMakeLists里做如下配置:

1
2
3
4
5
6
FIND_PACKAGE( OpenMP REQUIRED)
if(OPENMP_FOUND)
message("OPENMP FOUND")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
endif()

注意: OpenMP并不适合需要复杂的线程间同步和互斥的场合,这种情况下花的时间可能更长

gmapping使用OpenMP加速的语句: #pragma omp parallel for

  1. for循环的drawFromMotion之前

    1
    2
    3
    4
    5
    //#pragma omp parallel for
    for (unsigned int i = 0; i < m_particles.size(); i++) {

    m_particles[i].pose = m_motionModel.drawFromMotion(m_particles[i].pose, relPose,m_odoPose);
    }
  2. invalidateActiveArea之前

    1
    2
    3
    //#pragma omp parallel for
    for (unsigned int i = 0; i < m_particles.size(); i++) {
    m_matcher.invalidateActiveArea();
  3. scanMatch函数开头,实际就是对整个scanMatch并行化

    1
    2
    3
    inline void scanMatch(double * plainReading) {
    #pragma omp parallel for
    for (unsigned int i = 0; i < m_particles.size(); i++) {
  4. updateMap函数中

1
2
3
4
5
#pragma omp parallel for
for(int x=0; x < smap.getMapSizeX(); x++)
{
for(int y=0; y < smap.getMapSizeY(); y++)
{
  1. 重采样resample函数

对于保留下来的粒子进行更新,在并行化操作里面m_particles.push_back()会报错,因此需要把push_back()提出来,在外面的的for循环进行

1
2
3
4
5
6
7
8
9
10
11
#pragma omp parallel for
for(int i = 0; i<tmp_size;i++)
{
//对保留下来的粒子数据进行更新
//每个粒子的权重都设置为相同的值
temp[i].setWeight(0);
//为每个粒子更新running_scans
//增加了一帧激光数据 因此需要更新地图
m_matcher.registerScan(temp[i].map,temp[i].pose,plainReading);
//m_matcher.registerScan(temp[i].lowResolutionMap,temp[i].pose,plainReading);
}

为每个粒子更新地图时,同样可以并行化


(八) 源码分析6 scan match

先接上篇,不是第一帧激光数据的情况,上来是scanMatch函数,也就是GridSlamProcessor::scanMatch(在文件gridslamprocessor.hxx)。

主要部分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
inline void GridSlamProcessor::scanMatch(const double* plainReading)
{
// sample a new pose from each scan in the reference
/* 每个粒子都要进行 scan-match */
double sumScore=0;
int particle_number = m_particles.size();
//可以用openMP的方式来进行并行化,因此这里不能用迭代器,只能用下标的方式
//并行化之后,里面的循环会均匀分到各个不同的核里面
for (int i = 0; i < particle_number; i++)
{
OrientedPoint corrected;
double score, l, s;
/*进行scan-match 计算粒子的最优位姿,这是gmapping本来的做法*/
score = m_matcher.optimize(corrected, m_particles[i].map, m_particles[i].pose, plainReading);
/* 匹配成功则更新最优位姿 */
if (score>m_minimumScore)
{
m_particles[i].pose = corrected;
}
/* 扫描匹配不上,则使用里程计的数据,使用里程计数据不进行更新.
因为在进行扫描匹配之前 里程计已经更新过了*/
else
{
if (m_infoStream)
{
m_infoStream << "Scan Matching Failed, using odometry. Likelihood=" << l <<std::endl;
m_infoStream << "lp:" << m_lastPartPose.x << " " << m_lastPartPose.y << " "<< m_lastPartPose.theta <<std::endl;
m_infoStream << "op:" << m_odoPose.x << " " << m_odoPose.y << " "<< m_odoPose.theta;
}
}
// 粒子的最优位姿计算了之后,重新计算粒子的权重,相当于粒子滤波器中的观测步骤
/* 计算p(z|x,m),粒子的权重由粒子的似然来表示, 计算粒子的得分和权重(似然)
注意粒子的权重经过ScanMatch之后已经更新了
* 在论文中 粒子的权重不是用最优位姿的似然值来表示的
* 是用所有的似然值的和来表示的, s是得分 l是似然,也就是权重 */
m_matcher.likelihoodAndScore(s, l, m_particles[i].map, m_particles[i].pose, plainReading);
sumScore+=score;
m_particles[i].weight+= l;
m_particles[i].weightSum+= l;

//set up the selective copy of the active area
//by detaching the areas that will be updated
/*计算出来最优的位姿之后,进行地图的扩充 这里不会进行内存分配
*不进行内存分配的原因是这些粒子进行重采样之后有可能会消失掉,因此在后面进行重采样的时候统一进行内存分配。
*理论上来说,这里的操作是没有必要的,因为后面的重采样的时候还会进行一遍
*/
m_matcher.invalidateActiveArea();
m_matcher.computeActiveArea(m_particles[i].map, m_particles[i].pose, plainReading);
}
if (m_infoStream)
m_infoStream << "Average Scan Matching Score of all particles=" << sumScore/m_particles.size();
}

开始的都没意思,直接看ScanMatcher::optimize,在scanmatcher.cpp里,开头又是一个重要函数score,在scanmatcher.h里。它根据地图、激光数据、位姿迭代求解一个最优的新的位姿出来。

代码比较复杂,直接看整个scan函数(包括了optimize和score)的评分算法的理论:

对于栅格地图的环境地图模型,一般采用基于似然场模型(likelihood field range finder mode)的扫描匹配算法来进行激光雷达数据与地图的匹配。 该算法通过评估当前时刻位姿的激光雷达数据中每一个激光点与该时刻栅格地图的匹配程度,并将评估得到的每个激光点得分进行累加,得到当前时刻激光雷达数据与地图匹配的总分数,得分越高说明该时刻激光雷达信息与地图匹配程度越高,即机器人位姿越精确。由于该算法能够根据机器人当前时刻状态、 观测值和地图信息直接得到该时刻位姿的匹配程度,因而该算法常与基于粒子滤波的定位算法结合使用,用于选取各粒子中得分最高的位姿,并计算得分, 从而完成机器人位姿的确定,其原理如图:
基于似然场模型的扫描匹配示意图.png

黑色的格子代表前一时刻栅格地图中障碍物的位置, 圆形的点表示以里程计运动模型作为提议分布得到的机器人位姿估计为基础,将当前时刻激光雷达数据转换到栅格地图坐标系中的激光点的分布。把激光雷达的坐标转换到世界坐标系: 先旋转到机器人坐标系,然后再转换到世界坐标系。

该位姿下观测信息与地图匹配得分的具体算法步骤如下:对于当前激光雷达数据中任意一个激光点,设其端点在世界坐标系中坐标为 lp ,则其对应栅格地图中坐标为 phit ,对 phit周围八个栅格进行障碍物判定,并计算这些障碍物的平均坐标 pavr, 求出 pavr与 phit 的距离 dist 。 该激光点的评分可由下式表示
评分公式.png
ε为基于扫描匹配概率的激光方差,对该时刻所有激光点进行上式的计算,并将评分进行求和,分数越高说明该位姿的可信度越高。 对当前时刻所有粒子位姿进行扫描匹配运算, 并根据得分实现粒子的权重更新,最后求出粒子群的平均评分sumScore/m_particles.size()

如果设置的minimumScore参数过大,一些粒子的匹配会失败,因为要求太高了,gmapping会出现下列信息:
scan match失败.png
gmapping变成使用里程计进行位姿估计,这其实是不好的,因为从论文可知激光精度比里程计精确得多,但是注意粒子的权重计算还是调用likelihoodAndScore函数。

试验小强的scan match评分

直接修改gridslamprocessor.hxx的scanMatch函数。 小强的gmapping编译有问题,执行catkin_make --pkg gmapping之后,需要 cp /home/xiaoqiang/Documents/ros/devel/lib/gmapping/slam_gmapping /home/xiaoqiang/Documents/ros/src/gmapping/launch,否则roslaunch找不到节点文件

根据测试,激光评分多数在140~160,有时也会超过160.

for循环剩下的部分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
/*矫正成功则更新位姿*/
if (score>m_minimumScore)
{
m_particles[i].pose = corrected;
}
// 扫描匹配不上,则使用里程计的数据,使用里程计数据不进行更新,因为扫描匹配之前,里程计已经更新过了
else
{
//输出信息 这个在并行模式下可以会出现错位
if (m_infoStream)
{
m_infoStream << "Scan Matching Failed, using odometry. Likelihood=" << l <<endl;
m_infoStream << "lp:" << m_lastPartPose.x << " " << m_lastPartPose.y << " "<< m_lastPartPose.theta <<endl;
m_infoStream << "op:" << m_odoPose.x << " " << m_odoPose.y << " "<< m_odoPose.theta <<endl;
}
}
/* 粒子的最优位姿计算了之后,重新计算粒子的权重(相当于粒子滤波器中的观测步骤,
计算p(z|x,m)),粒子的权重由粒子的似然来表示。
* 计算粒子的得分和权重(似然) 注意粒子的权重经过ScanMatch之后已经更新了
* 在论文中,例子的权重不是用最有位姿的似然值来表示的。
* 是用所有的似然值的和来表示的。
*/
m_matcher.likelihoodAndScore(s, l, m_particles[i].map, m_particles[i].pose, plainReading);

sumScore+=score;
m_particles[i].weight+=l;
m_particles[i].weightSum+=l;
//set up the selective copy of the active area
//by detaching the areas that will be updated
/*计算出来最优的位姿之后,进行地图的扩充 这里不会进行内存分配
*不进行内存分配的原因是这些粒子进行重采样之后有可能会消失掉,因此在后面进行冲采样的时候统一进行内存分配。
*理论上来说,这里的操作是没有必要的,因为后面的重采样的时候还会进行一遍
*/
m_matcher.invalidateActiveArea();
m_matcher.computeActiveArea(m_particles[i].map, m_particles[i].pose, plainReading);

重采样

resample函数: 该函数在gridslamprocessor.hxx。首先是备份老的粒子的轨迹,即保留叶子的节点。然后是需要重采样还是不需要重采样,如果不需要重采样,则权值不变。只为轨迹创建一个新的节点,每个粒子更新地图。当有效值小于阈值的时候需要重采样,通过resampleIndexes提取到需要删除的粒子。删除粒子后,保留当前的粒子并在保存的粒子的节点里新增一个节点。删除要删除粒子的节点,保留的粒子进行数据更新,将每个粒子的设置同一个权重。最后更新一下地图。

resampleIndexes:该函数在particlefilter.h中,使用轮盘赌算法,决定哪些粒子会保留,保留的粒子会返回下标,里面的下标可能会重复,因为有些粒子会重复采样,而另外的一些粒子会消失掉。

首先计算总的权重,计算平均权重值(interval),根据权值进行采样,target是0-1分布随机选取的一数值,当总权重大于目标权重的,记录该粒子的索引,target在加上一个interval。如果某个粒子的权重比较大的话,就被采样多次。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
// 简化声明
int resampleIndexes(const std::vector<Particle>& particles, int nparticles) const
{
Numeric cweight=0;
/*计算总的权重*/
unsigned int n=0;
for (typename std::vector<Particle>::const_iterator it=particles.begin(); it!=particles.end(); ++it)
{
cweight+=(Numeric)*it;
n++;
}
if (nparticles>0)
n=nparticles;
//compute the interval
Numeric interval=cweight/n;

// drand48 返回服从均匀分布的·[0.0, 1.0) 之间的double随机数
Numeric target=interval*::drand48();
// 如果某个粒子的区间为4*interval。那么它至少被采样3次。
cweight=0; n=0;
std::vector<unsigned int> indexes(n);
unsigned int i=0;
for (typename std::vector<Particle>::const_iterator it=particles.begin(); it!=particles.end(); ++it, ++i)
{
cweight+=(Numeric)* it;
while(cweight>target)
{
indexes[n++]=i;
target+=interval;
}
}
return indexes;
}