Thursday 31 December 2015

How to debug garbage

I encounter several cases where the debugging infor gives me garbage like this

7f95d26c0000-7f95d26fd000 r-xp 00000000 08:05 5767281                    /lib/x86_64-linux-gnu/libpcre.so.3.13.1


It seems to be caused by wrong path, for example, there is dash (if there is space is ok). The path in windows is not always recognized by Ubantu. Correct the path format and this issue is solved. 

PCL pointer issue

I want to set a point cloud as a member of a class, and declare PointCloud::Ptr out_pcl_cloud; in the class definition. However, when I try to concatenate several point clouds together, it always gives me segmentation fault.

To solve this, declare a new pointer in the cpp file, PointCloud::Ptr cur_pcl_cloud (new PointCloud());, and use out_pcl_cloud = cur_pcl_cloud; to copy the point cloud data. Do not use *out_pcl_cloud = *cur_pcl_cloud; which appears to copy nothing.

To show all the pcl continuously, initialize the viewer by pcl::visualization::CloudViewer viewer( "viewer" );, and then viewer.showCloud( vo.out_pcl_cloud ); Do not use while( !viewer.wasStopped() ).

Tuesday 29 December 2015

ROS rviz resource not found

I use the launch file to start rviz and encounter the problem of "resource not found". According to here and here, this is because ros cannot find the package location. To solve the problem, go to the catkin workspace directory and source devel/setup.bash, and then go back to where the launch file is located, such as the sub folder of the package, and roslaunch.

Saturday 26 December 2015

opencv中用optical flow检测运动物体

这里描述了一个简单的方法可以在图像中检测运动物体。首先在两张图像中运行optical flow, 然后找到homography。之后用perspectiveTransform把prevPts投影到warpedPts,比较warpedPts和nextPts就能找到运动的物体。如果物体不动,warpedPts和nextPts很接近,不过如果有运动它们之间的差就会很大。

Friday 25 December 2015

一步步学习Visual odometry - Part2

有两种方式可以匹配两帧中的特征点,一是跟踪一个图像中的特征,适用于运动幅度较小的情况。二是在两张图像中分别检测特征点再进行匹配。为了去掉错误的匹配,可以用mutual consistency check和ratio test。不过ratio test需要根据经验设定阈值,所以不如用ransac来去掉outliers.

特征点的个数对结果有很大影响。一般特征点越多越好,而且要均匀分布在图像中,可以把图像分成网格,在每个格子中提取一定数目的特征点。对于640x480的图像,通常1000个特征点比较合适。

除了使用特征点,还可以用dense correspondence match。光流法可以跟踪全部的或者特定的像素,不过它假定运动很小,可能不适于VO。

用ransac去除outlier的时候,需要的内点越少所需要的iteration越少,时间也越少。所以尽可能的使得计算运动所需要的特征点数目比较小。不过如果特征点噪音很大,用更多的点会使得误差更加均匀分布,从而减小结果的误差。所以这时用8点算法可能好于5点算法。5点算法的code在opencv 3中有find essential mat函数,另外还有这里这里


一步步学习Visual odometry - Part1

在单目VO中,相机的运动和场景的3D结构都要从2D图像中获取。因为不知道绝对的距离,初始的两个相机之间的距离通常设为1.每当得到一个新的图像,相对的尺度和相机运动由已知的3D结构或者trifocal tensor得到。根据wiki,tensor是一个3x3x3的矩阵。trifocal tensor就是描述三帧之间关系的矩阵。

单目VO主要分为三类,基于特征点的方法,appearance based,和两者的混合。Nister第一次使用了VO这个词。跟之前方法不同的是,他们不是跟踪特征点,而是在每一帧里面单独检测Harris角点,然后进行匹配,这样就避免了跟踪过程中特征点的漂移。他们也不是用3D对3D的匹配,而是3D对2D的匹配,而且使用了5点算法和ransac去除outlier。还有一个是Milford提出的RatSLAM,用场景的中心做template tracking。

为了减小叠加相机姿态造成的轨迹漂移,通常在最近的m帧使用bundle adjustment优化,这个过程叫windowed BA。

在计算相机运动的时候,通过匹配特征点可以找到两帧之间的E矩阵。8点算法要求3D点不能都共面,而5点算法没有这个限制,它只需要相机标定。这里提到了Opencv中的findFundamentalMat使用了8点算法因此不是特别稳定,推荐使用5点算法。在该函数中有两个地方可以调整,一个是epipolar line的距离,一个是confidence level。根据这里,距离的设定是0.006 * maxVal。

在分解E矩阵的时候,可以得到四个可能的R,t。通过三角定位判断所有点是不是都在相机的前面,我们可以得到唯一解。这个解作为初始值再进行优化得到正确的R,t值。优化的过程就是减小反投影残差。

初始的两帧的运动确定之后,后续的尺度要以此为依据。前两帧的3D点对和对应的新的3D点对之间距离的比值就是尺度的比值。为了得到更准确的值,我们可以计算许多点的平均值或者中值。相应的,位移t也根据这个尺度调整。这里的特征点至少要在三帧中被跟踪和匹配。另外一种方法是用三帧之中特征点的trifocal constraint来计算尺度,不需要三角定位。




Thursday 24 December 2015

一步步学习Visual odometry - 综述

最近做视觉里程计visual odometry各种崩溃,只能使用一步步来的方法来梳理头绪。
Scaramuzza大神有关于VO很好的介绍,这里记录下学习笔记。
VO有几个基本假设。首先光照条件要充足,然后场景大部分要是静止的。图像中需要有足够的特征来提取运动。两帧之间的重叠要足够大。

VO最早用在火星车上面,后来被Nister发扬光大。它解决的问题是,一个相机在场景中运动,并记录下一系列图像,从这些图像中如何恢复出相机运动。把没两帧之间的运动相加,就可以得到一个完整的轨迹。我们也可以在过去的m帧上面运行iterative refinement来优化局部轨迹。

在VO中,运动估计是最重要的部分,它计算了相机从前一帧到当前帧的运动。我们可以用appearance based方法,直接用图像的像素值来匹配两帧图像。我们也可以用特征点来进行匹配。

如果两帧之间的匹配都是2D的点,至少需要5个点来恢复运动,点越多越好。通过这些对应点可以找到本质矩阵E,然后用SVD分解得到R和t。我们通过最小化反投影残差来找到真实的解。但是,我们并不知道在真实世界中相机移动的距离,只知道相机移动了一个单位。

如果我们既知道真实的3D点,又知道图像中对应的2D点,这就是PnP问题,至少需要3个点来解决,也是通过最小化反投影残差。在单目VO中,需要先从两帧之中通过三角定位得到3D点,然后跟第三帧的2D点匹配。

在三角定位中,如果相机移动的距离相比场景的深度很小,那么得到的3D点有很大的不确定性。为了避免这个情况,需要舍弃那些距离很近的帧,而是选取距离远的关键帧。选取的标准是3D点的不确定性很小。这一步是非常重要的,一定要在更新运动之前进行。

最后,在最近的m帧中要做bundle adjustment来优化轨迹。

在图像匹配时,会出现很多和运动不相符的outlier,我们必须把它们去掉。常用的方法是ransac。ransac不会给出确定的答案,所以不同时间的运行可能得到不同结果。当iteration增多时,得到的结果往往更稳定。最高效的方法是1点ransac,找到inlier后可以估算6DOF的运动。如果特征点的匹配有很大噪音,用更多的点的ransac会更准确。

为了优化VO的结果,通常需要上千个特征点,均匀分布在整个图像中。


Wednesday 23 December 2015

How to extract data from ROS bag

Extracting image from ros bag can use call back function directly and save the image. To save the data in text file, such as the camera matrix, rosbag play and then use the following command according to this

rostopic echo -p /camera/camera_info > data.txt

If you want to save a lot of data, such as the information from imu, see here

in the bag directory, run

rostopic echo -p -b airground_rig_s3_2013-03-18_21-38-48.bag /imu/data >imu.csv

and then you can open imu.csv in matlab.

Monday 21 December 2015

CV resources

I have found a good site here, which lists useful CV blogs and resources.

Sunday 20 December 2015

一步步学习SFM 五 Opencv multiple view geometry

之前做VO只能做到两帧之间对比的up to scale的版本,之后受到理论理解不深的限制很难再继续做优化之类的得到metric scale。看来还是需要好好先把理论搞清楚。正好看到了这本书,于是补一补基础知识。

分解F矩阵
两帧之中特征点和相机姿态的关系受到epipolar geometry的约束。通常假定第二个相机的姿态是相对第一个相机而言的,就是说第二个相机的旋转R和平移t都是相对第一个相机。两帧之中的特征点满足一个关系,叫做epipolar constrait,跟特征矩阵F有关。F又跟R和t有关。F的rank是2,determinant是0.也就是说相机的姿态可以从F中恢复,具体推导过程在这本书The essential and fundamental matrices这个章节。

假设一个3D点X在第一帧的投影是x1,在第二帧中的投影是x2,那么x1的位置使得x2只能在第二帧的一条线上,叫做epipolar line。F可以帮助在这条线上寻找x2的位置。所有的epipolar line都相交于一点,叫做epipole。

F是一个3x3的矩阵,不过尺度是不确定的,所以只需要有8个对应点就能计算F。如果特征点数目多于8个,就用least square找到误差最小的F。不过当8个点都在同一平面上的时候,这个方法不能使用。Opencv中的cv::findFundamentalMat函数还可以用7个点计算F,因为我们已经知道所有的epipolar line都交于一点,这个条件也限定了F。findFundamentalMat在用7点算法的时候可能返回3,因为这时可能有三个解。在8点算法,ransac, LMedS都只返回1或者0.检查findFundamentalMat的返回数很重要,因为有些时候点选取的不好,无法计算F。

这本书中提到了一个重要的问题,就是求解F的时候图像中对应点的选取非常关键。这些特征点需要在图像中均匀分布,并且包含不同深度的点。实际上可以用RANSAC来得到更准的F。首先用SURF寻找和描述特征点,再用knn matcher来找到两个最近的匹配点。如果最好匹配的距离比次好匹配的距离小很多,那么可以肯定最好匹配没有误差,否者说明最好和次好匹配很难区分,要去掉该特征点。接下来从图像A到B和图像B到A分别做匹配,只选取那些都是最好匹配的特征点。最后就是用RANSAC找到F和inlier。之前去掉错误匹配的过程是为了让ransca有个更好的初始条件,inlier更多就更容易得到准的F。

假如对应点都在同一个平面上,可以让它们的Z坐标都是0,那么3x4的projection matrix就变成了3x3的matrix,也就是homography。

当知道两个相机的姿态,可以用三角定位来确定一对对应特征点的3D坐标。如果已知3D点和图像中的投影点,也能算出相机姿态,这个问题叫做camera resectioning。

在已知相机矩阵的情况下,图像点的坐标根据相机矩阵K进行normalize之后,它们之间的关系跟essential matrix有关,这里用E矩阵来表示这是已知K的情况下的F。从E中可以计算出四个可能的相机姿态解,不过只有一个能让所有点都在两个相机前面。E和H的区别是,H假设了看到的点在同一个平面,所以可以把一个3D点定位到一个2D点。E没有平面假设,只能把一个3D点定位到另外一个图像的一条线上。E和F的区别是,E描述了同一个3D点在左右两个相机中的3D坐标的关系,而F描述了同一个2D点在左右两个相机中坐标的关系。

为了得到robust的解,还可以用ransac。在多帧的计算中相机姿态的误差会累积,这时用bundle adjustment优化3D点和相机姿态。

POSIT
这本书中提到了POSIT算法,用来计算已知物体大小的R,t运动。我们需要知道物体上4个不共面的点。算法首先假设这些点都在同一个深度,它们的大小变化仅仅跟相机距离物体远近有关。也就是说假设物体距离相机足够远,我们可以忽略物体的深度变化,所以物体上的点深度都一样。这个假设也叫weak perspective approximation.

从已知的点可以计算出一个并不是很准的pose,通过这个pose又可以反投影到4个2D点,然后根据它们和3D点的对应关系又可以更新pose,一直到pose收敛。

注意这里假设物体距离相机很远,所以物体上面点的深度变化可以忽略。如果这个假设不成立,那么pose要么是错误的,要么不会收敛。更多的共面的点也不会让算法更准,反而会导致错误的结果。更多的不共面的点则可以帮助减少误差。这个不是很理解,既然本来选取的点就不共面,表示深度不同,可是算法本身还假设深度相同,是不是有点矛盾?

在旧版的opencv中,cvPOSIT只能读取一个focal_length,因为这里假定像素的x和y方向是正方形。


Optical flow and velocity estimation

According to this, structure from motion solves two problems: correspondence for 2D motion in the image, and reconstruction for 3D motion of the camera.

Correspondence is often solved by optical flow. Let the pixel intensity be E(x,y,t), and suppose the brightness pattern is displaced by dx, dy over time dt, then the intensity of the displaced pixel should be the same as the original one: $E(x,y,t) = E(x+dx,y+dy,t+dt)$. In other words, the total derivative of E w.r.t time is zero: $\frac{\delta E}{\delta x}\frac{dx}{dt}+\frac{\delta E}{\delta y}\frac{dy}{dt}+\frac{\delta E}{\delta t} = 0$.

If we denote $E_x = \frac{\delta E}{\delta x}, E_y = \frac{\delta E}{\delta y}, E_t = \frac{\delta E}{\delta t}, v_x = \frac{dx}{dt}, v_y = \frac{dy}{dt}$, then the equation becomes $E_x v_x + E_y v_y + E_t = 0$, which is known as the brightness constancy equation (BCE). $E_x, E_y, E_t$ are measurable, while $v_x, v_y$ are the unknown flows. One BCE is insufficient to solve the flows. Therefore, to solve it we have to consider a small window of adjacent pixels. For instance, LK algorithm assumes the motion field at a given time is constant over a block of pixels. For the block of n pixels, there are n BCEs. This linear over-determined equation could be solved by least square estimation.

Once we have the flows, we can estimate the 3D motion. Homography describes the linear position relationship between feature points on two different planes. Suppose the camera takes images of the ground during flight, the visual features in two consecutive frames are related by a 3x3 matrix H, assuming the ground plane is flat. If R and T are the inter-frame rotation and translation, N is the unit normal vector of the ground plane in the camera frame, d is the altitude, then $H = R + \frac{1}{d}TN^T$, and thus $T=d(H-R)N$.

Based on this, when UAV is moving slowly, the rotation is negligible. Thus we can assume R is an identity matrix. As T is estimated, the velocity vector v in UAV body frame can be computed by $v = \frac{T}{\Delta t}$.

A simple linear KF could be used to estimate the UAV position. The measurement is the estimated velocity based on the homography, the linear states $x_k$ include the position and velocity of the UAV. The state x is updated using the prediction model: $x_k = F_k x_{k-1} + B_k u_k + w_k$. Essentially, it means $x_t = x_{t-1} + \dot{x_{t-1}} \delta t + \frac{1}{2}a \delta t^2 + N(0, Q_k)$ for position, and $\dot{x_t} = \dot{x_{t-1}} + a \delta t + N(0, Q_k)$. $u_k$ is the system input, i.e. the UAV acceleration. We assume the acceleration follows a zero mean normal distribution. $w_k$ is the process noise, which follows $N(0, Q_k)$.

The measurement model is $z_k = G_k x_k + v_k$, where $v_k$ is the measurement noise following $N(0, R_k)$.

OpenCV for secret agents - Motion amplifying camera

在007的电影中,Bond经常需要查看被隐藏的武器。通常他使用的x ray,红外线等仪器,但是这些仪器也很容易被其他间谍探测到。其实仅仅使用相机分析人的移动就能得到很多信息,比如某个人是不是在看一个方向,或者一条腿一直在动等等的肢体语言分析。我们可以把某种运动不断重复实现放大,就好比水面的波纹不断的扩散,这种方法叫做Eulerian video magnification.

Eulerian video magnification是基于Eulerian specification of the flow field。假设有一条河,Eulerian specificaiton描述了某个位置在某个时间的水流速度。这个速度在高山的泉水会快,在冬季的河口会慢,在河面会快,在河底会慢。我们可以用它来描述在某个位置,某个时间的物体运动。我们可以把不同时间和地点的运动叠加起来,达到强化运动的效果。

另外一种表达方式是Lagrangian specification,描述了某个物体在某段时间的位置,比如河里的一粒沙子肯呢过需要几年时间从山顶流到河口。它好比许多计算机视觉的程序,描述的是一个物体或者特征在一段时间内的移动。

FFT
我们通常用一段波形图来表示声音信号,波峰代表大声,波谷代表小声。在一段视频中,每一个像素也可以被看做是一个有起有落的波形。假设我们同时用相机和麦克风来记录一个节拍器,那么它发声的波形和指针在图像中运动的波形相一致。所以无论是声音还是图像都能用频率来描述。我们通过把图像转换到频域,选取特定的频率,再转换回来,保留想要看到的运动模式。

在原始的图像中运行FFT会很慢,而且这些频率对应在每个像素的运动,会有很多噪声,而且细节过多。因此我们把图像downsample,不过我们也希望保留对运动很重要的边缘。

可以用Gaussian filter进行downsampling,每个像素都是周围像素的加权平均,然后每隔一个像素去掉一个像素,这样每一帧的大小都是上一个图像的一半,就得到了Gaussian image pyramid。

为了保留边缘,我们用Laplacian image pyramid。在已经得到的Gaussian image pyramid中,我们选取第i+1层图像,通过复制像素值把它的大小扩大一倍,再使用Gaussian filter。我们把得到的图像和在Gaussian image pyramid第i层的图像相减,就得到了Laplacian image pyramid的第i层图像。所以Laplacian image pyramid中的图像就是一个模糊的,downsampled的图像和另外一个更加模糊的,经过两次downsample,一次upsample的图像的差值。

为啥这样就得到边缘了呢?边缘实际上是局部对比很大的特征,非边缘就是局部平滑的区域。如果我们把一个平滑区域变得模糊,它还是平滑的,所以和原来区域的差值为0.反之差值不为0,这样就能找到边缘。

我们不仅downsample FFT的输入,也upsample IFFT的输出,这样就可以得到和原始图像大小一样的图像。






Thursday 17 December 2015

Convert for_each to for

In c++11, you can use for_each as std::for_each(fts.begin(), fts.end(), [&](Feature* i){do something});, which is really convenient. However, this requires using -std=c++11 in cmake and it somehow conflicts with SSBA.

In order to convert for_each to for loop and avoiding to use -std=c++11, there are several things to do. Take the setExistingFeatures defined here for example,  we have to change const Features& fts to Features& fts first, otherwise an error of "conversion from const requested" will appear.

Next, for fts.begin(), fts.end(), we have to define an iterator, std::list<Feature*>::iterator it = fts.begin(); it != fts.end(); it++, which could be used in the for loop.

Finally note that it is a pointer to a structure, to access the members in the structure, we have to use (*it)->px[1].

OpenCV for secret agents - Rearview camera

这一章讲的是Bond的座驾如何自动检测到后面的追踪者。为了检测车辆后视镜里面的车头灯,可以假定它们都是白色的圆形,然后用blob detection,分为五步
1. 用二值化处理,把高于某个阈值的亮度设置为白色,反之为黑色。
2. 找到每个blob的轮廓。
3. 把相邻的blob融合在一起。
4. 计算blob的特征,比如中心点,半径等。
5. 去掉那些不符合条件的blob。

为了区分不同颜色的车灯,还需要转换到HSV颜色空间。
Hue指的是颜色在颜色轮上面的角度,红色是0, 绿色是120,蓝色是240。
Saturation指的是颜色与灰度值之间的距离,比如rgb的最大值减去最小值。
在Opencv中,为了能用一个byte存储hue的数值,hue被压缩到了0到180之间。想要更精确的数值,就要用0到360之间的浮点数。

在测量距离的时候,可以用 物体在图像中的大小/焦距 = 物体实际大小/距离。左边的单位可以是像素,右边是米。有几种方法可以测量距离
1. time of flight相机。这种相机发射一束光到物体上,然后测量反射光的强度。一些相机,比如Kinect,使用红外线光。但是这种相机会受到其他强烈光源的干扰,所以不适合检测车灯。
2. 双目相机由两个相机组成,它们之间的距离是已知的。通过图像每个点的disparity,就是像素距离,可以得到实际距离。不过通常计算量比较大。
3. Structure from motion, 只需要一个相机,不过需要知道它移动的距离。从不同地方的相机采集的图像中可以像双目相机那样算出距离。这里需要假设物体是静止的。因为007的座驾和后面的车灯都在移动,这种方法也不行。


Wednesday 16 December 2015

一步步学习SFM 三

之前写的VO在feature很丰富的时候,比如fountain系列,貌似可以使用了。但是在面对草地等特征不足的时候,或者帧数增加的时候,solvePnPRansac就会找不到足够的inlier。在这篇博客中,作者提到了关键帧提取的重要性。并不是每一帧都要加入地图中,而是选取距离较远而且匹配成功,运动估计也正确的关键帧加入到地图中。

在code中,关键帧的定义是vector< FRAME > keyframes;进行匹配后要判断NOT_MATCHED,TOO_FAR_AWAY, TOO_CLOSE。这里注意每一个新的图像都是跟上一个keyframe匹配,而不是跟上一个图像匹配。

为了检查为什么inlier不够,我们可以把反投影得到的点显示出来。
cv::line(reprojected,imgPoints[ppt],projected3D[ppt],cv::Scalar(0,0,255),1);表示在名叫reprojected的图中用红色的直线连接原始的点和反投影点。
cv::line(reprojected,imgPoints[inliers[ppt]],projected3D[inliers[ppt]],cv::Scalar(0,0,255),1);在图中显示inlier,也用红色的直线连接。
for(int ppt=0;ppt<imgPoints.size();ppt++) {
cv::circle(reprojected, imgPoints[ppt], 2, cv::Scalar(255,0,0), CV_FILLED);
cv::circle(reprojected, projected3D[ppt], 2, cv::Scalar(0,255,0), CV_FILLED);
}
用蓝色点表示原始点,用绿色点表示反投影点。
cv::circle(reprojected, imgPoints[inliers[ppt]], 2, cv::Scalar(255,255,0), CV_FILLED);用蓝色加上绿色表示inlier点。

在inlier非常少的时候,光流的匹配貌似没有问题,可能是Find2D3DCorrespondences出现了问题。不过即使在inlier非常少的时候,found 1187 3d-2d point correspondences, 说明对应点的数目并不少,不知道为什么solvePnPRansac无法计算R,t。

Find2D3DCorrespondences在之前的点云中寻找当前图像中也出现的匹配点,貌似没什么问题。可能是计算3D点的时候误差逐渐增大,导致一段时间之后2D和3D的对应虽然可以找到,但是计算出的运动完全错误。这个问题在点云中也可以看出来,多个图像之后点云更加分散,导致3D和2D匹配计算运动很不准。在TriangulatePointsBetweenViews中减小反投影点残差的阈值,可以得到更准的3D点云,得到的inlier也就更多。

无论误差大小,总是需要处理从新初始化的情况,所以要拼接点云。在Find2D3DCorrespondences之后是FindPoseEstimation得到R,t,然后TriangulatePointsBetweenViews。 这时候就要用到之前求得的相机姿态。

这里又发现一个问题,PruneMatchesBasedOnF和FindCameraMatrices都是在用GetFundamentalMat,只不过PruneMatchesBasedOnF多了一行matches_matrix[std::make_pair(working_view,older_view)] = FlipMatches(matches_matrix[std::make_pair(older_view,working_view)]);原始代码中两个函数都要运行。也许是8点算法需要先除掉outlier,然后用inlier再计算一遍F才能准。

因为暂时早不到问题的根源,所以先用两个图像匹配,利用kitti中提供的groundtruth作为尺度。奇怪的是运行图像0到7的时候,提示三角定位失败,但是单独运行图像6和7则是成功的。发现是TriangulatePoints中的P必须是identity matrix,不能是上一次相机的姿态。可能是因为上次的姿态计算有误差,造成三角定位的误差过大而出错。

除了运行的问题,它的速度也偏慢,只有2Hz.可以尝试使用更快的特征点检测,比如SVO中用的fast,具体安装参考这里

8算法不能处理所有点都在同一个平面的问题,具体见这里这里。为了更好的计算E,可以使用5点算法。算法的code可以在这里,和这里找到。

Tuesday 15 December 2015

OpenCV for secret agents - head gesture recogntion

本章主要讲解如何判断一个人在点头还是在摇头。人脸识别用的是opencv的cascadeclassifier,在识别成功后用goodFeaturesToTrack检测人脸的特征,然后用calcOpticalFlowPyrLK来追踪这些特征。如果特征点的中心上下移动,就认为是点头。如果是左右移动,就认为是摇头。

在跟踪特征的时候要用到光流,这里有三个附加说明:
1. 我们不区分相机和人脸的运动
2. 特征的颜色和亮度在图像中保持相似
3.相邻的特征有相似的运动

在opencv中,calcOpticalFlowPyrLK函数用了Lucas Kanade方法计算光流。这个方法基于每个特征点周围的$3\times 3$的图像块,然后在相邻图像中寻找能使得残差平方最小的匹配图像块。opencv中还用了图像金字塔来进行多尺度处理,这样就能处理大的或者小的运动。其他的函数,比如calcOpticalFlowSF或者calcOpticalFlowFarneback都是跟踪每一个点,所以calcOpticalFlowPyrLK的计算量更小,处理不同远近的人脸也更加稳定。

goodFeaturesToTrack选取了更容易被跟踪的点,抛弃了不易被跟踪的点,比如反射的阳光,树枝的交叉点,因为这些点会很快变化或者消失。相机的视角变换可以用warping来模拟,这样不同视角下都稳定出现的点就会被选取。good features to track是基于minimum eigenvalue corners和Harris corners的改进。更先进的特征点检测,比如FAST,ORB,SIFT,SURF,FREAK更广泛的选取特征点,计算量也更大,不适合用在基础的光流计算上。尤其是对于头部运动的识别来说,速度更快更重要,并不需要找到很多的特征点。

SVO论文笔记

SVO提出了一种基于向下看的相机和IMU的定位算法。论文开始先总结了运动估计的方法,主要有两种。

第一种是基于特征点,现在图像中选取角点等特征,然后在相邻图像中匹配,用epipolar集合来回复相机运动和场景。最后最小化反投影残差来优化相机姿态和场景。这种方法的成功在于有很多特征点检测和匹配的算法,甚至可以用在两个图像之间运动很大的情况。缺点是需要设定检测和匹配的阈值来剔除错误的匹配,而且很多检测方法往往不精确,导致运动的漂移需要平均很多特征点才能减小。

第二种是直接的方法,基于图像中的像素值。这里使用局部的梯度,包括大小和方向。当图像中特征点很少的时候,或者图像有模糊的时候,这种方法都好于特征点法。计算photometric error需要比反投影残差更多的计算量,不过不需要特征点的检测和匹配。

之前多数单目视觉里程计都是基于PTAM的,它是基于特征点检测的SLAM算法,通过跟踪和匹配上百个特征点实现。它通过并行计算运动估计和建图,使用高效的关键帧的BA,来实现实时性能。不过PTAM主要用在室内桌面的AR,如果在室外使用需要限制关键帧的个数。

直接的方法有DTAM,用最小化一个能量方程来计算每个关键帧里面的深度值。相机的姿态通过匹配整个深度图来计算。这种方法需要很大计算量,所以只能用GPU。

对比之前的方法,SVO主要的贡献在
1. semi-direct方法,比之前的方法更准确,快速。
2. 基于概率的构图方法,可以处理outlier。

SVO使用了两个并行的thread,一个用来估计相机姿态,另一个构建地图。

估计相机姿态时,第一步是最小化相邻两个图像中相同3D点投影的photometric error。反投影的点对应的2D点通过匹配图像块进行优化。最后同时优化姿态和反投影残差。

在论文的motion estimation部分,作者假设像素的残差成正态分布,variance是1.与之前方法不同,并不是很多的点都有深度,而是只有稀疏的特征点有深度。在每个特征点周围选取4x4的图像块,我们要找到能使得这些图像块残差最小的相机姿态。

为了减小漂移,相机的姿态应该和地图对应,而不是前一帧。所有能在当前相机中看到的已知3D点会被反投影到图像中,对每个反投影点,都找到它在关键帧中对应的图像块,然后在跟当前图像中的图像块对比。

在姿态和结构优化部分,作者对运动和3D点分别使用BA来优化,然后在将两者放在一起优化。跟只是两帧之间对比的sparse image alignment比较,加入当前图像和关键帧以及地图的对比能得到更准确的结果。

在地图构建中,我们根据当前图像和相机的姿态来估计每个特征点的深度。构建地图的时候,一个基于概率的深度滤波器会对每个2D特征点初始化,来寻找对应3D点。开始的时候每个深度的不确定性都非常大,之后每个图像中都根据贝叶斯定理优化深度。当一个深度的不确定性足够小的时候,一个新的3D点就被加入地图,用来做运动估计。

这部分可以参考这篇论文,大概是讲如何从多个图像中估计深度。在图像I中,在光轴上一点x(Z),距离相机深度为Z,它在I中的图像块w对应另外一个图像I'中的图像块w',可以用NCC来计算w和w'的相似度。通过观察作者发现用寻找NCC最大值的方法计算Z,要么得到的深度很准,要么就随机分布在Zmin和Zmax的区间内,也就是可以用一个gaussian+uniform的概率分布来描述。给定一个初始的先验概率,接下来就可以用bayesian inference来计算最大化后验概率的深度。

在SVO的论文里,depth filter的大意貌似是这样的。计算了某个点u在图像Ir的深度rP之后,在图像Ik中可以找到u的对应点u',如果把u'移动一个像素,就能得到一个新的u对应的3D点rP+,对深度估计的误差就是rP和rP+之间的距离。为什么rP的variance是把u'移动一个像素之后得到的呢?具体要参考这里。更多关于为何使用depth filter参考这里

仔细看mapping部分,每个depth filter都对应一个关键帧。初始化的时候每个点的不确定性很大,初始值就是关键帧中深度的平均值。接下来的每个图像中,我们在epipolar line上面寻找对应的图像块,然后通过三角定位计算深度。跟这篇论文中不同的是,为了应对很大的深度变化,SVO使用了inverse depth。

算法实现
在最开始的两帧,SVO使用了平面假设来估计一个homography,然后用三角定位初始化两帧里面点的深度,这个过程和PTAM相似。为了应对大幅度的运动,在sparse image alignment用了coarse to fine方法。

如果当前帧距离所有关键帧的距离超过场景平均深度的12%,那么当前帧就被选为关键帧。为了加快效率,SVO中的关键帧数目是一定的,所以每加入一个关键帧,都会把最老的一个关键帧去掉。

在mapping thread中,图像被分成了30x30的小格子。每个Shi-Tomasi score最高的FAST角点会被加入depth filter,除非那个格子里面已经有了2D和3D对应的特征点存在。这样就可以得到平均分布的特征点。

推荐另外一篇非常好的SVO论文博客

运行SVO
从github上面下载并安装好svo之后,可以运行test pipeline来看结果。在catkin_ws里面catkin_make, 生成的.exe文件在catkin_ws/devel/lib/svo里面。作者给的airground dataset的pose可以用这里提到的方法存在txt文件里面。








Monday 14 December 2015

Visual odometry - OpenMVO

I have found several great sources on VO, namely this, this, and this. Since I would like to work on monocular VO, I use this first. In this blog, I will write about its implementation.

The original code is written for windows, and I want to port it to Ubantu in ROS catkin workspace. So I do not follow the file structure but use the one for catkin. Whenever a third party library is needed, I include it in the cmake.

At first, we have to detect features using FAST. We need to define a class called frame to process the images. Several third party libraries are needed, similar to SVO. To use Eigen, I copy the CMakeModules folder in SVO to my src folder, which contains FindEigen.cmake. Regarding Sophus, I encounter an error saying that Sophus.cmake cannot be found. This is wired because in SVO it finds the library without any problem. After debugging, it turns out that I should use FIND_PACKAGE(Sophus REQUIRED) in cmake instead of FIND_PACKAGE(sophus REQUIRED). The name should be exactly the same and the first letter should be CAPITAL.

For the constructor in the frame class, : key_pts_(5) means to set key_pts_ to 5, as mentioned here.  In Frame.cpp, it turns out that we could not use #include "feature.h", which leads to re-definition of frame class.

We also have to add

IF(CMAKE_COMPILER_IS_GNUCC)
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x")
ELSE()
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
ENDIF()

as in SVO cmake to be able to use lambda expression, eg. 
std::for_each(key_pts_.begin(), key_pts_.end(), [&](Feature* ftr){ ftr = NULL; });

Since I am using Ubantu, I could not find _aligned_malloc, which seems to be only available in windows VS. So I delete the SSE2 related part in image_utils.h. 

Another issue often encountered is the re-declaration problem. For instance, in the point3d.h, we need to add
#ifndef POINT3D_H_
#define POINT3D_H_
and #endif
to avoid this.

In frame.h, when it says Vector3d does not name a type, we have to include abstract_camera.h, and  using namespace Eigen; solves this problem.

When I run the initial_position, it gives core dumped error, and it seems c2f is not executed. It also turns out that px_ref_it->x - px_cur_it->x, px_ref_it->y - px_cur_it->y are all 0. This is because I use two identical images for the frame 1 and 2. But even if I use different images, it still gives me segmentation fault. The reason is that I forget to initialize cam_(cam) in frame.cpp

Take note that vector is in namespace std, not cv.

Next problem is that scale is nan. The implementation seems to be problematic, and I probably should implement my own VO step by step.

In the FAST keypoint detection, the image is divided into several cells, and only one keypoint in every cell with the highest shi-tomasi score is selected, as shown in the following code
if (score > corners.at(k).score)
corners.at(k) = Corner(xy.x*scale, xy.y*scale, score, L, 0.0f);
There is an angle in the Corner defined in feature.h, which should hold the gradient information of the keypoint, but it does not seem to be in the code.  

OpenCV for secret agents - cat recogntion

The cat recognition chapter in the book reminds me of the recent 007 movie, Specter. The interesting things to note in this chapter is that the training cat images have to be rotated to upright position. The function getRotationMatrix2D is used to define the rotation and warpAffine is used to rotate the image. Obviously some black regions will be introduced in the image, and thus we can fill them by a specific color. The author uses 50% gray, as it does not tend to introduce bias into the image histogram equalization.


Haar cascades and LBPH

Opencv for secret agents是一本非常棒的opencv指导用书。有一章讲的是machine learning。在一张图像中我们可以看到有特点的形状或者模式,而这些模式往往只在特定尺度下出现。为了生成Haar cascade,图像的不同区域会经历尺度变化,这样我们就能每次只关注一些像素。像这样对图像取样叫做window。我们可以用window中的一些区域的像素减去另外一些区域的像素值来判断是否包含一个灰暗的部分和一个明亮的部分,从而可以对比当前window和一些常见形状,比如边缘和角点。如果一个window和某个特定形状很相似,我们就可以选择它作为一个feature。对于同一物体,通常我们能够在相似的地方和尺度找到相似的feature。
并非所有的feature都同样显著。如果一个feature在物体中很常见,我们就把它放到positive training set里面。反之,不常见的就放到negative training set里面。根据feature描述物体的程度,我们可以赋给它们不同的stage.多个stage就组成一个cascade。为了判断待检测目标是不是已知样本,每个stage都必须通过。

Thursday 3 December 2015

一步步学习SFM 二

上一篇看了理论,这里主要看code。

在读取所有图像后, 用MultiCameraPnP初始化。如果没有给相机矩阵,在MultiCameraDistance里面会自动初始化一个。

OnlyMatchFeatures会用OFFeatureMatcher把所有已经读取的图像之间的光流计算出来并进行匹配,然后把结果存储在matches_matrix里面。

然后PruneMatchesBasedOnF会计算F,然后把不符合的匹配去掉。GetBaseLineTriangulation可以在所有图像中找到两张homography匹配最好的,然后计算R,t。


AdjustCurrentBundle使用SBA把点云和相机轨迹进行优化。其中有个参数是K = Count2DMeasurements(pointcloud);就是同一个点在几个图像中都被观察到了。如果没有检查新得到的点是否在已经有的点云中,不知道还能不能使用BA。


如果只计算当前frame与上一张之间的光流,在计算F的时候使用的也是这两个frame的匹配。FindCameraMatrices会输出tmp_pcloud表示outCloud。

这里遇到一个问题,就是TriangulatePointsBetweenViews总是找不到原来点云中的匹配点,导致所有的点都被认为是新的。在onlymatchfeature中,每一张新的frame都要跟原来所有的frame做匹配,否则TriangulatePointsBetweenViews找不到对应点。TriangulatePointsBetweenViews是以older view和working view为基准,所以不能只让当前frame和前一frame进行匹配。

还有一个问题是GetRGBForPointCloud总是报错,BUG: point id:1235490744 should not exist for img #2 which has only 9280,函数是在AdjustCurrentBundle中调用的GetRGBForPointCloud(pointcloud_beforeBA,pointCloudRGB_beforeBA);

可能是int pt_idx = _pcloud[i].imgpt_for_img[good_view]超过边界的问题。经过检查,_pcloud[i].imgpt_for_img.size()是2,表示imgpt_for_img里面只有working view和older view的匹配。后来发现PruneMatchesBasedOnF只处理了当前frame和前一个frame,加入for之后_pcloud[i].imgpt_for_img.size()就是3了。但是那个BUG还是会出现。

经过长时间debug,终于发现是TriangulatePointsBetweenViews函数中初始化new_triangulated[j].imgpt_for_img = std::vector<int>(2,-1);的问题,imgpt_for_img的size可能比实际用到的小。

should be pcloud[pt3d].imgpt_for_img[1] == submatches[ii].queryIdx instead of pcloud[pt3d].imgpt_for_img[0] == submatches[ii].queryIdx

在使用code的时候,如果是fountain dataset效果还可以,换成其他图像可能会计算不出有效F,需要把FindCameraMatrices函数中if(fabsf(determinant(E)) > 1e-07)阈值调高。

在连续运行的时候,会出现not enough inliers to consider a good pose,是因为在solvePnPRansac中inlier不够。根据这里的分析,可能是没有关键帧的提取。之后我会分析如何改进。