ORB-SLAM内的卡方检验

Reference:

1. 卡方检验(Chi-square test/Chi-Square Goodness-of-Fit Test)
2. 卡方检验详解分析与实例

1. 概念

• 卡方值: ${\chi }^2$ 值表示观察值与理论值之间的偏离程度。计算这种偏离程度的基本思路如下：

  1. 设 $O$ 代表某个类别的观察频数，$E$ 代表基于某个假设 $H_0$ 计算出的期望频数，$O$ 与 $E$ 之差称为残差；
  2. 残差可以表示某一个类别观察值和理论值的偏离程度，但如果将残差简单相加以表示各类别观察频数与期望频数的差别，则有一定的不足之处。因为残差有正有负，相加后会彼此抵消，综合仍然为 $0$，为此可以将残差平方后求和；
  3. 另一方面，残差大小是一个相对概念，相对于期望频数为 $10$ 时，期望频数为 $20$ 的残差非常大，但相对于期望频数为 $1000$ 时 $20$ 的残差就很小了。考虑到这一点，人们又将残差平方除以期望频数再求和，以估计观察频数与期望频数的差别。

  进行上述操作之后，就得到了常用的 ${\chi }^2$ 统计量，其计算公式为：
  $${\chi }^2=\sum \frac{(O-E{)}^2}{E}=\sum _{i=1}^k\frac{(O_i-E_i{)}^2}{E_i}=\sum _{i=1}^k\frac{(O_i-n{p}_i{)}^2}{n{p}_i}(\mathrm{i}=1,2,3,...,\mathrm{k})$$其中，$O_i$ 为 $i$ 水平的观察频数，$E_i$ 为 $i$ 水平的期望频数，$n$ 为总频数，$p_i$ 为 $i$ 水平的期望概率，$i$ 水平的期望频数 $E_i$ 等于总频数 $n\times i$ 水平的期望概率 $p_i$，$k$ 为单元格数。当 $n$ 比加大时，${\chi }^2$ 统计量近似服从 $k-1$(计算 $E_i$ 时用到的参数个数)个自由度的卡方分布。

  由卡方的计算公式可知，当观察频数与期望频数完全一致时，${\chi }^2$ 值为 $0$；观察频数与期望频数越接近，两者之间的差异越小，${\chi }^2$ 值越小；反之，观察频数与期望频数差别越大，两者之间的差异越大，${\chi }^2$ 值越大。换言之，大的 ${\chi }^2$ 值表明观察频数远离期望频数，即表明远离假设。小的 ${\chi }^2$ 值表明观察频数接近期望频数，接近假设。因此，${\chi }^2$ 是观察频数与期望频数之间距离的一种度量指标，也是假设成立与否的度量指标。如果 ${\chi }^2$ 值小，研究者就倾向于不拒绝 $H_0$；如果 ${\chi }^2$ 值大，就倾向于拒绝 $H_0$。至于 ${\chi }^2$ 在每个具体研究中究竟要达到什么程度才能拒绝 $H_0$，则要借助于卡方分布求出所对应的 $P$ 来确定，这就引出了显著性水平的概念。

• 显著性水平: 显著性水平是估计总体参数落在某一区间内可能犯错的概率，用 $\alpha$ 表示，即当原假设为正确时，却把它滤掉了的概率或风险，通常取 $\alpha =0.05$ 或 $\alpha =0.01$。

2. 卡方检验的基本思想

卡方检验是以 ${\chi }^2$ 分布为基础的一种常用假设检验方法，它的无效建设 $H_0$ 是：观察频数与期望频数没有差别。

该检验的基本思想是：首先假设 $H_0$ 成立，基于此前提计算出 ${\chi }^2$ 的值，它表示观察值与理论值之间的偏离程度。根据 ${\chi }^2$ 分布及自由度可以确定在 $H_0$ 假设成立的情况下获得当前统计量及更极端情况的概率 $P$。如果 $P$ 值很小，说明观察值与理论值偏离程度太大，应当拒绝无效假设，表示比较资料之间有显著差异；否则就不能拒绝无效假设，尚不能认为样本所代表的实际情况和理论假设有差别。

根据自由度和显著性水平查询检验统计量临界值，其中纵轴为自由度，横轴为显著性水平阈值。比如检查单应阵的函数，点到点的重投影误差自由度为 $2$，在显著性水平为 $0.05$ 时通过查表得阈值为 $5.99$(需要注意的是，卡方值为 ${\chi }^2$，而不是单独的 $\chi$)：

3. 卡方检测示例

为了验证肺癌与吸烟的关系，我们得到如下数据：

是否肺病患者	吸烟	不吸烟	合计	吸烟比例
是	158	169	327	48%
否	82	311	393	20%
合计	240	480	720	33%

首先假设吸烟与肺癌两者之间没有关系(这里应该不用先假设有没关系，越大就是越有关系，越小就是越没有关系，根据数值判断相关性)，我们计算期望值(因为上面吸烟的比例为33%，因此在吸烟与肺癌没有关系的时候，肺癌患者吸烟与不吸烟的比例应该是33%，没有得肺癌的吸烟与不吸烟的比例也应该是33%)：

是否肺病患者	吸烟	不吸烟	合计	吸烟比例
是	109	218	327	33%
否	131	262	393	33%
合计	240	480	720	33%

带入卡方计算公式：
$${\chi }^2=\frac{(158-109{)}^2}{109}+\frac{(169-218{)}^2}{218}+\frac{(82-131{)}^2}{131}+\frac{(311-262{)}^2}{262}=60.53$$自由度的计算方法，可以简单抽象成：(行数-1)*(列数-1)，所以这里的自由度为 $1$。

通过查表可得自由度为 $1$ 时，显著性水平为 $0.05$，当卡方值小于 $3.84$ 时，可以接受原假设，即变量之间没有相关性。卡方值越小，不相关的概率越大。现在卡方值远大于 $3.84$，说明两者不相关的概率很小，即抽烟与肺癌有关。

4. ORB-SLAM2中卡方检测剔除外点的策略

就特征点法的视觉SLAM而言，一般会计算重投影误差。具体而言，记 $\mathbf{u}$ 为特征点的 $2D$ 位置，$\bar{\mathbf{u}}$ 为由地图点投影到图像上的 $2D$ 位置。重投影误差为：
$$\mathbf{e}=\mathbf{u}-\bar{\mathbf{u}}$$重投影误差服从高斯分布：
$$\mathbf{e}\sim \mathcal{N}(0,\mathbf{\Sigma })$$其中，协方差 $\mathbf{\Sigma }$ 一般根据特征点提取的金字塔层级决定。具体的，记提取ORB特征时，图像金字塔的每层缩小尺度为 $s$(ORB-SLAM中为1.2)。在ORB-SLAM中假设第 $0$ 层的标准差为 $p$ 个pixel(ORB-SLAM中设为了1个pixel)；那么，一个在金字塔第 $n$ 层提取的特征的重投影误差的协方差为：
$$\mathbf{\Sigma }=(s^n\times p{)}^2\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$式 $\mathbf{e}=\mathbf{u}-\bar{\mathbf{u}}$ 中的误差是一个 $2$ 维向量，这里阈值不好设置，就把它变成一个标量，计算向量的内积 $r$(向量元素的平方和)。但是，这又有一个问题，不同金字塔层的特征点都用同一个阈值，是不是不合理呢？于是，在计算内积的时候，利用协方差进行加权(协方差表达了不确定度)。金字塔层级越高特征点提取精度越低，协方差 $\mathbf{\Sigma }$ 越大，那么就有了：
$$r={\mathbf{e}}^T{\mathbf{\Sigma }}^{-1}\mathbf{e}$$利用协方差加权，起到了归一化的作用，具体如上式，可以变为：
$$r=({\mathbf{\Sigma }}^{-\frac{1}{2}}\mathbf{e}{)}^T({\mathbf{\Sigma }}^{-\frac{1}{2}}\mathbf{e})$$而：
$$({\mathbf{\Sigma }}^{-\frac{1}{2}}\mathbf{e})\sim \mathcal{N}(0,\mathbf{I})$$为多维标准正态分布，也就是说不同金字塔层提取的特征，计算的重投影误差都被归一化了，或者说去量纲化了，那么，我们只用一个阈值就可以了。

说白了就是计算协方差加权的重投影误差和卡方作为阈值间的大小，若误差超过了阈值，则该被当成外点：$\text{score}=e^T\Sigma e>{\chi }^2?\text{outlier}:\text{inlier}$

4.1 示例，卡方检验计算置信度得分: CheckFundamental()、CheckHomography()

根据重投影误差构造统计量 ${\chi }^2$，其值越大，观察结果和期望结果之间的差别越显著，某次计算越可能用到了外点。

说白了其实就是 乘上信息矩阵的平方误差：

/**
 * @brief 对给定的homography matrix打分,需要使用到卡方检验的知识
 * 
 * @param[in] H21                       从参考帧到当前帧的单应矩阵
 * @param[in] H12                       从当前帧到参考帧的单应矩阵
 * @param[in] vbMatchesInliers          匹配好的特征点对的Inliers标记
 * @param[in] sigma                     方差，默认为1
 * @return float                        返回得分
 */
float Initializer::CheckHomography(
    const cv::Mat &H21,                 //从参考帧到当前帧的单应矩阵
    const cv::Mat &H12,                 //从当前帧到参考帧的单应矩阵
    vector<bool> &vbMatchesInliers,     //匹配好的特征点对的Inliers标记
    float sigma)                        //估计误差
{
    // 说明：在已知n维观测数据误差服从N(0，sigma）的高斯分布时
    // 其误差加权最小二乘结果为  sum_error = SUM(e(i)^T * Q^(-1) * e(i))
    // 其中：e(i) = [e_x,e_y,...]^T, Q维观测数据协方差矩阵，即sigma * sigma组成的协方差矩阵
    // 误差加权最小二次结果越小，说明观测数据精度越高
    // 那么，score = SUM((th - e(i)^T * Q^(-1) * e(i)))的分数就越高
    // 算法目标： 检查单应变换矩阵
    // 检查方式：通过H矩阵，进行参考帧和当前帧之间的双向投影，并计算起加权最小二乘投影误差

    // 算法流程
    // input: 单应性矩阵 H21, H12, 匹配点集 mvKeys1
    //    do:
    //        for p1(i), p2(i) in mvKeys:
    //           error_i1 = ||p2(i) - H21 * p1(i)||2
    //           error_i2 = ||p1(i) - H12 * p2(i)||2
    //           
    //           w1 = 1 / sigma / sigma
    //           w2 = 1 / sigma / sigma
    // 
    //           if error1 < th
    //              score +=   th - error_i1 * w1
    //           if error2 < th
    //              score +=   th - error_i2 * w2
    // 
    //           if error_1i > th or error_2i > th
    //              p1(i), p2(i) are inner points
    //              vbMatchesInliers(i) = true
    //           else 
    //              p1(i), p2(i) are outliers
    //              vbMatchesInliers(i) = false
    //           end
    //        end
    //   output: score, inliers

	// 特点匹配个数
    const int N = mvMatches12.size();

	// Step 1 获取从参考帧到当前帧的单应矩阵的各个元素
    const float h11 = H21.at<float>(0,0);
    const float h12 = H21.at<float>(0,1);
    const float h13 = H21.at<float>(0,2);
    const float h21 = H21.at<float>(1,0);
    const float h22 = H21.at<float>(1,1);
    const float h23 = H21.at<float>(1,2);
    const float h31 = H21.at<float>(2,0);
    const float h32 = H21.at<float>(2,1);
    const float h33 = H21.at<float>(2,2);

	// 获取从当前帧到参考帧的单应矩阵的各个元素
    const float h11inv = H12.at<float>(0,0);
    const float h12inv = H12.at<float>(0,1);
    const float h13inv = H12.at<float>(0,2);
    const float h21inv = H12.at<float>(1,0);
    const float h22inv = H12.at<float>(1,1);
    const float h23inv = H12.at<float>(1,2);
    const float h31inv = H12.at<float>(2,0);
    const float h32inv = H12.at<float>(2,1);
    const float h33inv = H12.at<float>(2,2);

	// 给特征点对的Inliers标记预分配空间
    vbMatchesInliers.resize(N);

	// 初始化score值
    float score = 0;

    // 基于卡方检验计算出的阈值（假设测量有一个像素的偏差）
	// 自由度为2的卡方分布，显著性水平为0.05，对应的临界阈值
    const float th = 5.991;

    //信息矩阵，方差平方的倒数
    const float invSigmaSquare = 1.0/(sigma * sigma);

    // Step 2 通过H矩阵，进行参考帧和当前帧之间的双向投影，并计算起加权重投影误差
    // H21 表示从img1 到 img2的变换矩阵
    // H12 表示从img2 到 img1的变换矩阵 
    for(int i = 0; i < N; i++)
    {
		// 一开始都默认为Inlier
        bool bIn = true;

		// Step 2.1 提取参考帧和当前帧之间的特征匹配点对
        const cv::KeyPoint &kp1 = mvKeys1[mvMatches12[i].first];
        const cv::KeyPoint &kp2 = mvKeys2[mvMatches12[i].second];
        const float u1 = kp1.pt.x;
        const float v1 = kp1.pt.y;
        const float u2 = kp2.pt.x;
        const float v2 = kp2.pt.y;

        // Step 2.2 计算 img2 到 img1 的重投影误差
        // x1 = H12*x2
        // 将图像2中的特征点通过单应变换投影到图像1中
        // |u1|   |h11inv h12inv h13inv||u2|   |u2in1|
        // |v1| = |h21inv h22inv h23inv||v2| = |v2in1| * w2in1inv
        // |1 |   |h31inv h32inv h33inv||1 |   |  1  |
		// 计算投影归一化坐标
        const float w2in1inv = 1.0/(h31inv * u2 + h32inv * v2 + h33inv);
        const float u2in1 = (h11inv * u2 + h12inv * v2 + h13inv) * w2in1inv;
        const float v2in1 = (h21inv * u2 + h22inv * v2 + h23inv) * w2in1inv;
   
        // 计算重投影误差 = ||p1(i) - H12 * p2(i)||2
        const float squareDist1 = (u1 - u2in1) * (u1 - u2in1) + (v1 - v2in1) * (v1 - v2in1);
        const float chiSquare1 = squareDist1 * invSigmaSquare;

        // Step 2.3 用阈值标记离群点，内点的话累加得分
        if(chiSquare1>th)
            bIn = false;    
        else
            // 误差越大，得分越低
            score += th - chiSquare1;

        // 计算从img1 到 img2 的投影变换误差
        // x1in2 = H21*x1
        // 将图像2中的特征点通过单应变换投影到图像1中
        // |u2|   |h11 h12 h13||u1|   |u1in2|
        // |v2| = |h21 h22 h23||v1| = |v1in2| * w1in2inv
        // |1 |   |h31 h32 h33||1 |   |  1  |
		// 计算投影归一化坐标
        const float w1in2inv = 1.0/(h31*u1+h32*v1+h33);
        const float u1in2 = (h11*u1+h12*v1+h13)*w1in2inv;
        const float v1in2 = (h21*u1+h22*v1+h23)*w1in2inv;

        // 计算重投影误差 
        const float squareDist2 = (u2-u1in2)*(u2-u1in2)+(v2-v1in2)*(v2-v1in2);
        const float chiSquare2 = squareDist2*invSigmaSquare;
 
        // 用阈值标记离群点，内点的话累加得分
        if(chiSquare2>th)
            bIn = false;
        else
            score += th - chiSquare2;   

        // Step 2.4 如果从img2 到 img1 和 从img1 到img2的重投影误差均满足要求，则说明是Inlier point
        if(bIn)
            vbMatchesInliers[i]=true;
        else
            vbMatchesInliers[i]=false;
    }
    return score;
}