[SLAM] GMapping SLAM源码阅读(草稿)

目前可以从很多地方得到RBPF的代码,主要看的是Cyrill Stachniss的代码,据此进行理解。

Author:Giorgio Grisetti; Cyrill Stachniss  http://openslam.org/ 

https://github.com/Allopart/rbpf-gmapping   和文献[1]上结合的比较好,方法都可以找到对应的原理。

https://github.com/MRPT/mrpt MRPT中可以采用多种扫描匹配的方式,可以通过配置文件进行配置。


 

双线程和程序的基本执行流程

GMapping采用GridSlamProcessorThread后台线程进行数据的读取和运算,在gsp_thread.h和相应的实现文件gsp_thread.cpp中可以看到初始化,参数配置,扫描数据读取等方法。

最关键的流程是在后台线程调用的方法void * GridSlamProcessorThread::fastslamthread(GridSlamProcessorThread* gpt)

而这个方法中最主要的处理扫描数据的过程在428行,bool processed=gpt->processScan(*rr);同时可以看到GMapping支持在线和离线两种模式。

注意到struct GridSlamProcessorThread : public GridSlamProcessor ,这个后台线程执行类GridSlamProcessorThread 继承自GridSlamProcessor,所以执行的是GridSlamProcessor的processScan方法。

  1 //后台线程处理扫描数据的方法
  2 bool GridSlamProcessor::processScan(const RangeReading & reading, int adaptParticles)
  3 {
  4     /**retireve the position from the reading, and compute the odometry*/
  5     OrientedPoint relPose=reading.getPose();
  6     if (!m_count)
  7     {
  8         m_lastPartPose=m_odoPose=relPose;
  9     }
 10     
 11     //write the state of the reading and update all the particles using the motion model
 12     for (ParticleVector::iterator it=m_particles.begin(); it!=m_particles.end(); it++)
 13     {
 14         OrientedPoint& pose(it->pose);
 15         pose=m_motionModel.drawFromMotion(it->pose, relPose, m_odoPose);//运动模型,更新t时刻的粒子
 16     }
 17 
 18     // update the output file
 19     if (m_outputStream.is_open())
 20     {
 21         m_outputStream << setiosflags(ios::fixed) << setprecision(6);
 22         m_outputStream << "ODOM ";
 23         m_outputStream << setiosflags(ios::fixed) << setprecision(3) << m_odoPose.x << " " << m_odoPose.y << " ";
 24         m_outputStream << setiosflags(ios::fixed) << setprecision(6) << m_odoPose.theta << " ";
 25         m_outputStream << reading.getTime();
 26         m_outputStream << endl;
 27     }
 28     if (m_outputStream.is_open())
 29     {
 30         m_outputStream << setiosflags(ios::fixed) << setprecision(6);
 31         m_outputStream << "ODO_UPDATE "<< m_particles.size() << " ";
 32         for (ParticleVector::iterator it=m_particles.begin(); it!=m_particles.end(); it++)
 33         {
 34             OrientedPoint& pose(it->pose);
 35             m_outputStream << setiosflags(ios::fixed) << setprecision(3) << pose.x << " " << pose.y << " ";
 36             m_outputStream << setiosflags(ios::fixed) << setprecision(6) << pose.theta << " " << it-> weight << " ";
 37         }
 38         m_outputStream << reading.getTime();
 39         m_outputStream << endl;
 40     }
 41     
 42     //invoke the callback
 43     onOdometryUpdate();
 44     
 45     // accumulate the robot translation and rotation
 46     OrientedPoint move=relPose-m_odoPose;
 47     move.theta=atan2(sin(move.theta), cos(move.theta));
 48     m_linearDistance+=sqrt(move*move);
 49     m_angularDistance+=fabs(move.theta);
 50     
 51     // if the robot jumps throw a warning
 52     if (m_linearDistance>m_distanceThresholdCheck)
 53     {
 54         cerr << "***********************************************************************" << endl;
 55         cerr << "********** Error: m_distanceThresholdCheck overridden!!!! *************" << endl;
 56         cerr << "m_distanceThresholdCheck=" << m_distanceThresholdCheck << endl;
 57         cerr << "Old Odometry Pose= " << m_odoPose.x << " " << m_odoPose.y 
 58         << " " <<m_odoPose.theta << endl;
 59         cerr << "New Odometry Pose (reported from observation)= " << relPose.x << " " << relPose.y 
 60         << " " <<relPose.theta << endl;
 61         cerr << "***********************************************************************" << endl;
 62         cerr << "** The Odometry has a big jump here. This is probably a bug in the   **" << endl;
 63         cerr << "** odometry/laser input. We continue now, but the result is probably **" << endl;
 64         cerr << "** crap or can lead to a core dump since the map doesn't fit.... C&G **" << endl;
 65         cerr << "***********************************************************************" << endl;
 66     }
 67     
 68     m_odoPose=relPose;
 69     
 70     bool processed=false;
 71 
 72     // process a scan only if the robot has traveled a given distance
 73     if (! m_count || m_linearDistance>m_linearThresholdDistance || m_angularDistance>m_angularThresholdDistance)
 74     {     
 75         if (m_outputStream.is_open())
 76         {
 77             m_outputStream << setiosflags(ios::fixed) << setprecision(6);
 78             m_outputStream << "FRAME " <<  m_readingCount;
 79             m_outputStream << " " << m_linearDistance;
 80             m_outputStream << " " << m_angularDistance << endl;
 81         }
 82       
 83         if (m_infoStream)
 84             m_infoStream << "update frame " <<  m_readingCount << endl
 85                 << "update ld=" << m_linearDistance << " ad=" << m_angularDistance << endl;
 86       
 87       
 88         cerr << "Laser Pose= " << reading.getPose().x << " " << reading.getPose().y 
 89         << " " << reading.getPose().theta << endl;
 90       
 91       
 92         //this is for converting the reading in a scan-matcher feedable form
 93         assert(reading.size()==m_beams);
 94         double * plainReading = new double[m_beams];
 95         for(unsigned int i=0; i<m_beams; i++)
 96         {
 97             plainReading[i]=reading[i];
 98         }
 99         m_infoStream << "m_count " << m_count << endl;
100         if (m_count>0)
101         {
102             scanMatch(plainReading);//扫描匹配
103             if (m_outputStream.is_open())
104             {
105                 m_outputStream << "LASER_READING "<< reading.size() << " ";
106                 m_outputStream << setiosflags(ios::fixed) << setprecision(2);
107                 for (RangeReading::const_iterator b=reading.begin(); b!=reading.end(); b++)
108                 {
109                     m_outputStream << *b << " ";
110                 }
111                 OrientedPoint p=reading.getPose();
112                 m_outputStream << setiosflags(ios::fixed) << setprecision(6);
113                 m_outputStream << p.x << " " << p.y << " " << p.theta << " " << reading.getTime()<< endl;
114                 m_outputStream << "SM_UPDATE "<< m_particles.size() << " ";
115                 for (ParticleVector::const_iterator it=m_particles.begin(); it!=m_particles.end(); it++)
116                 {
117                     const OrientedPoint& pose=it->pose;
118                     m_outputStream << setiosflags(ios::fixed) << setprecision(3) <<  pose.x << " " << pose.y << " ";
119                     m_outputStream << setiosflags(ios::fixed) << setprecision(6) <<  pose.theta << " " << it-> weight << " ";
120                 }
121                 m_outputStream << endl;
122             }
123             onScanmatchUpdate();
124             updateTreeWeights(false);//更新权重
125                 
126             if (m_infoStream)
127             {
128                 m_infoStream << "neff= " << m_neff  << endl;
129             }
130             if (m_outputStream.is_open())
131             {
132                 m_outputStream << setiosflags(ios::fixed) << setprecision(6);
133                 m_outputStream << "NEFF " << m_neff << endl;
134             }
135             resample(plainReading, adaptParticles);//重采样
136         } 
137         else 
138         {
139             m_infoStream << "Registering First Scan"<< endl;
140             for (ParticleVector::iterator it=m_particles.begin(); it!=m_particles.end(); it++)
141             {    
142                 m_matcher.invalidateActiveArea();
143                 m_matcher.computeActiveArea(it->map, it->pose, plainReading);//计算有效区域
144                 m_matcher.registerScan(it->map, it->pose, plainReading);
145       
146                 // cyr: not needed anymore, particles refer to the root in the beginning!
147                 TNode* node=new    TNode(it->pose, 0., it->node,  0);
148                 node->reading=0;
149                 it->node=node;
150       
151             }
152         }
153         //cerr  << "Tree: normalizing, resetting and propagating weights at the end..." ;
154         updateTreeWeights(false);//再次更新权重
155         //cerr  << ".done!" <<endl;
156       
157         delete [] plainReading;
158         m_lastPartPose=m_odoPose; //update the past pose for the next iteration
159         m_linearDistance=0;
160         m_angularDistance=0;
161         m_count++;
162         processed=true;
163       
164         //keep ready for the next step
165         for (ParticleVector::iterator it=m_particles.begin(); it!=m_particles.end(); it++)
166         {
167             it->previousPose=it->pose;
168         }
169       
170     }
171     if (m_outputStream.is_open())
172         m_outputStream << flush;
173     m_readingCount++;
174     return processed;
175 }
GridSlamProcessor::processScan

可以依次看到下面介绍的1-7部分的内容。

 

1.运动模型

t时刻粒子的位姿最初由运动模型进行更新。在初始值的基础上增加高斯采样的noisypoint,参考MotionModel::drawFromMotion()方法。原理参考文献[1]5.4.1节,sample_motion_model_odometry算法。

p是粒子的t-1时刻的位姿,pnew是当前t时刻的里程计读数,pold是t-1时刻的里程计读数。参考博客:小豆包的学习之旅:里程计运动模型

 1 OrientedPoint MotionModel::drawFromMotion(const OrientedPoint& p, const OrientedPoint& pnew, const OrientedPoint& pold) const{
 2     double sxy=0.3*srr;
 3     OrientedPoint delta=absoluteDifference(pnew, pold);
 4     OrientedPoint noisypoint(delta);
 5     noisypoint.x+=sampleGaussian(srr*fabs(delta.x)+str*fabs(delta.theta)+sxy*fabs(delta.y));
 6     noisypoint.y+=sampleGaussian(srr*fabs(delta.y)+str*fabs(delta.theta)+sxy*fabs(delta.x));
 7     noisypoint.theta+=sampleGaussian(stt*fabs(delta.theta)+srt*sqrt(delta.x*delta.x+delta.y*delta.y));
 8     noisypoint.theta=fmod(noisypoint.theta, 2*M_PI);
 9     if (noisypoint.theta>M_PI)
10         noisypoint.theta-=2*M_PI;
11     return absoluteSum(p,noisypoint);
12 }

 

2.扫描匹配

扫描匹配获取最优的采样粒子。GMapping默认采用30个采样粒子。

 1 /**Just scan match every single particle.
 2 If the scan matching fails, the particle gets a default likelihood.*/
 3 inline void GridSlamProcessor::scanMatch(const double* plainReading){
 4   // sample a new pose from each scan in the reference
 5   double sumScore=0;
 6   for (ParticleVector::iterator it=m_particles.begin(); it!=m_particles.end(); it++){
 7     OrientedPoint corrected;
 8     double score, l, s;
 9     score=m_matcher.optimize(corrected, it->map, it->pose, plainReading);
10     //    it->pose=corrected;
11     if (score>m_minimumScore){
12       it->pose=corrected;
13     } else {
14     if (m_infoStream){
15       m_infoStream << "Scan Matching Failed, using odometry. Likelihood=" << l <<std::endl;
16       m_infoStream << "lp:" << m_lastPartPose.x << " "  << m_lastPartPose.y << " "<< m_lastPartPose.theta <<std::endl;
17       m_infoStream << "op:" << m_odoPose.x << " " << m_odoPose.y << " "<< m_odoPose.theta <<std::endl;
18     }
19     }
20 
21     m_matcher.likelihoodAndScore(s, l, it->map, it->pose, plainReading);
22     sumScore+=score;
23     it->weight+=l;
24     it->weightSum+=l;
25 
26     //set up the selective copy of the active area
27     //by detaching the areas that will be updated
28     m_matcher.invalidateActiveArea();
29     m_matcher.computeActiveArea(it->map, it->pose, plainReading);
30   }
31   if (m_infoStream)
32     m_infoStream << "Average Scan Matching Score=" << sumScore/m_particles.size() << std::endl;    
33 }
GridSlamProcessor::scanMatch

注意ScanMatcher::score()函数的原理是likehood_field_range_finder_model方法,参考《概率机器人》手稿P143页。

ScanMatcher::optimize()方法获得了一个最优的粒子,基本流程是按照预先设定的步长不断移动粒子的位置,根据提议分布计算s,得到score最小的那个作为粒子的新的位姿。

 1 //此处的方法是likehood_field_range_finder_model方法,参考《概率机器人》手稿P143
 2 inline double ScanMatcher::score(const ScanMatcherMap& map, const OrientedPoint& p, const double* readings) const{
 3     double s=0;
 4     const double * angle=m_laserAngles+m_initialBeamsSkip;
 5     OrientedPoint lp=p;
 6     lp.x+=cos(p.theta)*m_laserPose.x-sin(p.theta)*m_laserPose.y;
 7     lp.y+=sin(p.theta)*m_laserPose.x+cos(p.theta)*m_laserPose.y;
 8     lp.theta+=m_laserPose.theta;
 9     unsigned int skip=0;
10     double freeDelta=map.getDelta()*m_freeCellRatio;
11     for (const double* r=readings+m_initialBeamsSkip; r<readings+m_laserBeams; r++, angle++){
12         skip++;
13         skip=skip>m_likelihoodSkip?0:skip;
14         if (*r>m_usableRange) continue;
15         if (skip) continue;
16         Point phit=lp;
17         phit.x+=*r*cos(lp.theta+*angle);
18         phit.y+=*r*sin(lp.theta+*angle);
19         IntPoint iphit=map.world2map(phit);
20         Point pfree=lp;
21         pfree.x+=(*r-map.getDelta()*freeDelta)*cos(lp.theta+*angle);
22         pfree.y+=(*r-map.getDelta()*freeDelta)*sin(lp.theta+*angle);
23          pfree=pfree-phit;
24         IntPoint ipfree=map.world2map(pfree);
25         bool found=false;
26         Point bestMu(0.,0.);
27         for (int xx=-m_kernelSize; xx<=m_kernelSize; xx++)
28         for (int yy=-m_kernelSize; yy<=m_kernelSize; yy++){
29             IntPoint pr=iphit+IntPoint(xx,yy);
30             IntPoint pf=pr+ipfree;
31             //AccessibilityState s=map.storage().cellState(pr);
32             //if (s&Inside && s&Allocated){
33                 const PointAccumulator& cell=map.cell(pr);
34                 const PointAccumulator& fcell=map.cell(pf);
35                 if (((double)   )> m_fullnessThreshold && ((double)fcell )<m_fullnessThreshold){
36                     Point mu=phit-cell.mean();
37                     if (!found){
38                         bestMu=mu;
39                         found=true;
40                     }else
41                         bestMu=(mu*mu)<(bestMu*bestMu)?mu:bestMu;
42                 }
43             //}
44         }
45         if (found)
46             s+=exp(-1./m_gaussianSigma*bestMu*bestMu);//高斯提议分布
47     }
48     return s;
49 }
View Code

ScanMatcher::likelihoodAndScore()和ScanMatcher::score()方法基本一致,但是是将新获得的粒子位姿进行计算q,为后续的权重计算做了准备。

ScanMatcher::optimize()方法——粒子的运动+score()中激光扫描观测数据。

其他的扫描匹配方法也是可以使用的:比如ICP算法(rbpf-gmapping的使用的是ICP方法,先更新了旋转角度,然后采用提议分布再次优化粒子)、Cross Correlation、线特征等等。

 

3.提议分布 (Proposal distribution)

注意是混合了运动模型和观测的提议分布,将扫描观测值整合到了提议分布中[2](公式9)。因此均值和方差的计算与单纯使用运动模型作为提议分布的有所区别。提议分布的作用是获得一个最优的粒子,同时用来计算权重,这个体现在ScanMatcher::likelihoodAndScore()和ScanMatcher::score()方法中,score方法中采用的是服从0均值的高斯分布近似提议分布(文献[1] III.B节)。关于为什么采用混合的提议分布,L(i)区域小的原理文献[1]中III.A节有具体介绍。

rbpf-gmapping的使用的是运动模型作为提议分布。

 

4.权重计算

在重采样之前进行了一次权重计算updateTreeWeights(false);

重采样之后又进行了一次。代码在gridslamprocessor_tree.cpp文件中。

1 void  GridSlamProcessor::updateTreeWeights(bool weightsAlreadyNormalized)
2 {
3   if (!weightsAlreadyNormalized)
4   {
5     normalize();
6   }
7   resetTree();
8   propagateWeights();
9 }

 

5.重采样

原理:粒子集对目标分布的近似越差,则权重的方差越大。

因此用Neff作为权重值离差的量度。

  1 //判断并进行重采样
  2 inline bool GridSlamProcessor::resample(const double* plainReading, int adaptSize, const RangeReading* )
  3 {
  4     bool hasResampled = false;
  5     TNodeVector oldGeneration;
  6     for (unsigned int i=0; i<m_particles.size(); i++)
  7     {
  8         oldGeneration.push_back(m_particles[i].node);
  9     }
 10     if (m_neff<m_resampleThreshold*m_particles.size())
 11     {
 12         if (m_infoStream)
 13             m_infoStream  << "*************RESAMPLE***************" << std::endl;
 14     
 15         uniform_resampler<double, double> resampler;
 16         m_indexes=resampler.resampleIndexes(m_weights, adaptSize);//执行重采样
 17     
 18         if (m_outputStream.is_open())
 19         {
 20             m_outputStream << "RESAMPLE "<< m_indexes.size() << " ";
 21             for (std::vector<unsigned int>::const_iterator it=m_indexes.begin(); it!=m_indexes.end(); it++)
 22             {
 23                 m_outputStream << *it <<  " ";
 24             }
 25             m_outputStream << std::endl;
 26         } 
 27         onResampleUpdate();
 28         //BEGIN: BUILDING TREE
 29         ParticleVector temp;//重采样的粒子集合临时变量
 30         unsigned int j=0;
 31         std::vector<unsigned int> deletedParticles;          //this is for deleteing the particles which have been resampled away.
 32     
 33         //cerr << "Existing Nodes:" ;
 34         for (unsigned int i=0; i<m_indexes.size(); i++)
 35         {
 36             //cerr << " " << m_indexes[i];
 37             while(j<m_indexes[i])
 38             {
 39                 deletedParticles.push_back(j);
 40                 j++;
 41             }
 42             if (j==m_indexes[i])
 43             j++;
 44             Particle & p=m_particles[m_indexes[i]];
 45             TNode* node=0;
 46             TNode* oldNode=oldGeneration[m_indexes[i]];
 47             //cerr << i << "->" << m_indexes[i] << "B("<<oldNode->childs <<") ";
 48             node=new TNode(p.pose, 0, oldNode, 0);
 49             node->reading=0;
 50             //cerr << "A("<<node->parent->childs <<") " <<endl;
 51       
 52             temp.push_back(p);
 53             temp.back().node=node;
 54             temp.back().previousIndex=m_indexes[i];
 55         }
 56         while(j<m_indexes.size())
 57         {
 58             deletedParticles.push_back(j);
 59             j++;
 60         }
 61         //cerr << endl;
 62         std::cerr <<  "Deleting Nodes:";
 63         //删除粒子
 64         for (unsigned int i=0; i<deletedParticles.size(); i++)
 65         {
 66             std::cerr <<" " << deletedParticles[i];
 67             delete m_particles[deletedParticles[i]].node;
 68             m_particles[deletedParticles[i]].node=0;
 69         }
 70         std::cerr  << " Done" <<std::endl;
 71     
 72         //END: BUILDING TREE
 73         std::cerr << "Deleting old particles..." ;
 74         m_particles.clear();
 75         std::cerr << "Done" << std::endl;
 76         std::cerr << "Copying Particles and  Registering  scans...";
 77         for (ParticleVector::iterator it=temp.begin(); it!=temp.end(); it++)
 78         {
 79             it->setWeight(0);
 80             m_matcher.invalidateActiveArea();
 81             m_matcher.registerScan(it->map, it->pose, plainReading);
 82             m_particles.push_back(*it);
 83         }
 84         std::cerr  << " Done" <<std::endl;
 85         hasResampled = true;
 86     }
 87     else 
 88     {
 89         int index=0;
 90         std::cerr << "Registering Scans:";
 91         TNodeVector::iterator node_it=oldGeneration.begin();
 92         for (ParticleVector::iterator it=m_particles.begin(); it!=m_particles.end(); it++)
 93         {
 94             //create a new node in the particle tree and add it to the old tree
 95             //BEGIN: BUILDING TREE  
 96             TNode* node=0;
 97             node=new TNode(it->pose, 0.0, *node_it, 0);
 98       
 99             node->reading=0;
100             it->node=node;
101 
102             //END: BUILDING TREE
103             m_matcher.invalidateActiveArea();
104             m_matcher.registerScan(it->map, it->pose, plainReading);//注册到全局地图中
105             it->previousIndex=index;
106             index++;
107             node_it++;   
108         }
109         std::cerr  << "Done" <<std::endl; 
110     }
111     //END: BUILDING TREE
112     return hasResampled;
113 }
GridSlamProcessor::resample

 

6.占用概率栅格地图

 此处的方法感觉有点奇怪,在resample方法中执行ScanMatcher::registerScan()方法,计算占用概率栅格地图。采样两种方式,采用信息熵的方式和文献[1] 9.2节的计算方法不一样。

 1 double ScanMatcher::registerScan(ScanMatcherMap& map, const OrientedPoint& p, const double* readings)
 2 {
 3     if (!m_activeAreaComputed)
 4         computeActiveArea(map, p, readings);
 5         
 6     //this operation replicates the cells that will be changed in the registration operation
 7     map.storage().allocActiveArea();
 8     
 9     OrientedPoint lp=p;
10     lp.x+=cos(p.theta)*m_laserPose.x-sin(p.theta)*m_laserPose.y;
11     lp.y+=sin(p.theta)*m_laserPose.x+cos(p.theta)*m_laserPose.y;
12     lp.theta+=m_laserPose.theta;//激光器中心点
13     IntPoint p0=map.world2map(lp);//转到地图坐标?
14     
15     
16     const double * angle=m_laserAngles+m_initialBeamsSkip;
17     double esum=0;
18     for (const double* r=readings+m_initialBeamsSkip; r<readings+m_laserBeams; r++, angle++)//每一条扫描线
19         if (m_generateMap)
20         {
21             double d=*r;
22             if (d>m_laserMaxRange)
23                 continue;
24             if (d>m_usableRange)
25                 d=m_usableRange;
26             Point phit=lp+Point(d*cos(lp.theta+*angle),d*sin(lp.theta+*angle));//扫描到的点
27             IntPoint p1=map.world2map(phit);//转到地图坐标?
28             IntPoint linePoints[20000] ;
29             GridLineTraversalLine line;
30             line.points=linePoints;
31             GridLineTraversal::gridLine(p0, p1, &line);//计算扫描线占据的栅格?
32             for (int i=0; i<line.num_points-1; i++)
33             {
34                 PointAccumulator& cell=map.cell(line.points[i]);
35                 double e=-cell.entropy();
36                 cell.update(false, Point(0,0));
37                 e+=cell.entropy();
38                 esum+=e;
39             }
40             if (d<m_usableRange)
41             {
42                 double e=-map.cell(p1).entropy();
43                 map.cell(p1).update(true, phit);
44                 e+=map.cell(p1).entropy();
45                 esum+=e;
46             }
47         } 
48         else
49         {
50             if (*r>m_laserMaxRange||*r>m_usableRange) continue;
51             Point phit=lp;
52             phit.x+=*r*cos(lp.theta+*angle);
53             phit.y+=*r*sin(lp.theta+*angle);
54             IntPoint p1=map.world2map(phit);
55             assert(p1.x>=0 && p1.y>=0);
56             map.cell(p1).update(true,phit);
57         }
58     //cout  << "informationGain=" << -esum << endl;
59     return esum;
60 }
ScanMatcher::registerScan

rbpf-gmapping的使用的是文献[1] 9.2节的计算方法,在Occupancy_Grid_Mapping.m文件中实现,同时调用Inverse_Range_Sensor_Model方法。

gridlineTraversal实现了beam转成栅格的线。对每一束激光束,设start为该激光束起点,end为激光束端点(障碍物位置),使用Bresenham划线算法确定激光束经过的格网。

小豆包的学习之旅:占用概率栅格地图和cost-map

 

7.计算有效区域

第一次扫描,count==0时,如果激光观测数据超出了范围,更新栅格地图的范围。同时确定有效区域。

  1 void ScanMatcher::computeActiveArea(ScanMatcherMap& map, const OrientedPoint& p, const double* readings){
  2     if (m_activeAreaComputed)
  3         return;
  4     OrientedPoint lp=p;
  5     lp.x+=cos(p.theta)*m_laserPose.x-sin(p.theta)*m_laserPose.y;
  6     lp.y+=sin(p.theta)*m_laserPose.x+cos(p.theta)*m_laserPose.y;
  7     lp.theta+=m_laserPose.theta;
  8     IntPoint p0=map.world2map(lp);
  9     
 10     Point min(map.map2world(0,0));
 11     Point max(map.map2world(map.getMapSizeX()-1,map.getMapSizeY()-1));
 12            
 13     if (lp.x<min.x) min.x=lp.x;
 14     if (lp.y<min.y) min.y=lp.y;
 15     if (lp.x>max.x) max.x=lp.x;
 16     if (lp.y>max.y) max.y=lp.y;
 17     
 18     /*determine the size of the area*/
 19     const double * angle=m_laserAngles+m_initialBeamsSkip;
 20     for (const double* r=readings+m_initialBeamsSkip; r<readings+m_laserBeams; r++, angle++){
 21         if (*r>m_laserMaxRange) continue;
 22         double d=*r>m_usableRange?m_usableRange:*r;
 23         Point phit=lp;
 24         phit.x+=d*cos(lp.theta+*angle);
 25         phit.y+=d*sin(lp.theta+*angle);
 26         if (phit.x<min.x) min.x=phit.x;
 27         if (phit.y<min.y) min.y=phit.y;
 28         if (phit.x>max.x) max.x=phit.x;
 29         if (phit.y>max.y) max.y=phit.y;
 30     }
 31     //min=min-Point(map.getDelta(),map.getDelta());
 32     //max=max+Point(map.getDelta(),map.getDelta());
 33     
 34     if ( !map.isInside(min)    || !map.isInside(max)){
 35         Point lmin(map.map2world(0,0));
 36         Point lmax(map.map2world(map.getMapSizeX()-1,map.getMapSizeY()-1));
 37         //cerr << "CURRENT MAP " << lmin.x << " " << lmin.y << " " << lmax.x << " " << lmax.y << endl;
 38         //cerr << "BOUNDARY OVERRIDE " << min.x << " " << min.y << " " << max.x << " " << max.y << endl;
 39         min.x=( min.x >= lmin.x )? lmin.x: min.x-m_enlargeStep;
 40         max.x=( max.x <= lmax.x )? lmax.x: max.x+m_enlargeStep;
 41         min.y=( min.y >= lmin.y )? lmin.y: min.y-m_enlargeStep;
 42         max.y=( max.y <= lmax.y )? lmax.y: max.y+m_enlargeStep;
 43         map.resize(min.x, min.y, max.x, max.y);
 44         //cerr << "RESIZE " << min.x << " " << min.y << " " << max.x << " " << max.y << endl;
 45     }
 46     
 47     HierarchicalArray2D<PointAccumulator>::PointSet activeArea;
 48     /*allocate the active area*/
 49     angle=m_laserAngles+m_initialBeamsSkip;
 50     for (const double* r=readings+m_initialBeamsSkip; r<readings+m_laserBeams; r++, angle++)
 51         if (m_generateMap)
 52         {
 53             double d=*r;
 54             if (d>m_laserMaxRange)
 55                 continue;
 56             if (d>m_usableRange)
 57                 d=m_usableRange;
 58             Point phit=lp+Point(d*cos(lp.theta+*angle),d*sin(lp.theta+*angle));
 59             IntPoint p0=map.world2map(lp);
 60             IntPoint p1=map.world2map(phit);
 61             
 62             IntPoint linePoints[20000] ;
 63             GridLineTraversalLine line;
 64             line.points=linePoints;
 65             GridLineTraversal::gridLine(p0, p1, &line);
 66             for (int i=0; i<line.num_points-1; i++)
 67             {
 68                 assert(map.isInside(linePoints[i]));
 69                 activeArea.insert(map.storage().patchIndexes(linePoints[i]));
 70                 assert(linePoints[i].x>=0 && linePoints[i].y>=0);
 71             }
 72             if (d<m_usableRange){
 73                 IntPoint cp=map.storage().patchIndexes(p1);
 74                 assert(cp.x>=0 && cp.y>=0);
 75                 activeArea.insert(cp);
 76             }
 77         } 
 78         else 
 79         {
 80             if (*r>m_laserMaxRange||*r>m_usableRange) continue;
 81             Point phit=lp;
 82             phit.x+=*r*cos(lp.theta+*angle);
 83             phit.y+=*r*sin(lp.theta+*angle);
 84             IntPoint p1=map.world2map(phit);
 85             assert(p1.x>=0 && p1.y>=0);
 86             IntPoint cp=map.storage().patchIndexes(p1);
 87             assert(cp.x>=0 && cp.y>=0);
 88             activeArea.insert(cp);
 89         }
 90     
 91     //this allocates the unallocated cells in the active area of the map
 92     //cout << "activeArea::size() " << activeArea.size() << endl;
 93 /*    
 94     cerr << "ActiveArea=";
 95     for (HierarchicalArray2D<PointAccumulator>::PointSet::const_iterator it=activeArea.begin(); it!= activeArea.end(); it++){
 96         cerr << "(" << it->x <<"," << it->y << ") ";
 97     }
 98     cerr << endl;
 99 */        
100     map.storage().setActiveArea(activeArea, true);
101     m_activeAreaComputed=true;
102 }
computeActiveArea

每次扫描匹配获取t时刻的最优粒子后会计算有效区域。

重采样之后,调用ScanMatcher::registerScan() 方法,也会重新计算有效区域。

 

参考文献:

[1]Sebastian Thrun et al. "Probabilistic Robotics(手稿)."

[2]Grisetti, G. and C. Stachniss "Improved Techniques for Grid Mapping with Rao-Blackwellized Particle Filters."

posted @ 2016-07-01 19:36  太一吾鱼水  阅读(29988)  评论(8编辑  收藏  举报