改善深度相机获取的深度图质量（稳定性）的trick

用kinect或者realsense获取的深度图或者点云图，由于传感器获取深度信息的不稳定性，经常出现一张图里面，有些深度信息丢失，尤其是物体边缘，因此怎么进行预处理得到较好的质量的数据很关键。

（自用备忘，若侵则删。）

1.时序数据的统计方法

多取N个时刻的数据（因为相机和物体保持不变），得到N张深度图，理论上这些数据应该相同，但由于传感器的不稳定性，有小部分不一致，因此可通过N张深度图，其每个pixel的值可进行累计求均值来改善，大致算法如下（来自别人的硕士论文“基于RealSense的散乱零件三维目标识别”）

2.双边滤波和联合双边滤波

为保持边沿特征，在滤波时要考虑核范围内的值，不能只用空间阈值。

双边滤波：

 

联合双边滤波

联合双边滤波相对于双边滤波来说，最大的特点就是引入了一幅引导影像。去这幅引导滤波灰度值来计算色域权重，引导图像要跟原图像一样大小，如果用了一幅不一样的图像做引导，需要进行上采样或者下采样，上采样的结果可以对深度图填补一些空洞的作用。

3.引导滤波

可参考，csdn的一篇文章

4.利用卡尔曼滤波算法

之前没了解卡尔曼滤波，看到了这个应用，就去了解了一下，看了B站的视频和CSDN的一篇博文（卡尔曼详解_上善若水-CSDN博客_卡尔曼），

基本公式：

　X(k)=AX(k-1)+BU(k)+W(k)       系统方程   （1）

　Z(k)=HX(k)+V(k)                         观测方程 （2）  

　X(k)是k时刻的系统状态，U(k)是k时刻对系统的控制量。A和B是系统参数，对于多模型系统，他们为矩阵。Z(k)是k时刻的测量值，H是测量系统的参数，对于多测量系统，H为矩阵。W(k)和V(k)分别表示过程和测量的噪声。他们被假设成高斯白噪声，他们的协方差分别是Q，R，在实际应用中。由于系统中一般不太有控制量，所以B这个参数一般为0，也就是没有U(K)（如果是自动驾驶，有加速度之类的，这些就是控制量）。           

注：虽然在表达式中包括了系统噪声w和转测噪声v，但在计算系统状态转移和系统观测量的时候并不考虑这两个噪声的影响。两个噪声的影响只是在卡尔曼滤波器离散迭代算法过程中使用到了两个噪声的协方差矩阵Q和R。分别用于计算系统状态估计误差的协方差矩阵P和卡尔曼滤波器增益K的大小。                           

自己的大致总结：这是世界任何数据或多或少都会有误差的存在，没有绝对的，我们看到的数据（头脑反馈估计数据-模型估计），设备测量出来的数据（传感器测出的数据），这些误差都是正态分布的。对于深度图来说，上一帧的优化后数据（每个pixel的值）为X（k-1），在下一帧预测模型（观察者大脑反馈）预测的值为X（k），我们的大脑反馈认为这一帧跟上一帧是一样的，所以A为1，B为0），这一帧的传感器测出来的是Z（k），就是得到的这一帧的深度图数据（每个pixel的值），这时候由于各种误差，X（k）和Z（k）并不是一样的（比如在传感器对物体边缘丢失，就是我们看到的黑点，本来就是有的，或者反过来），那怎么去权衡，（到底那个值是真是不知道的，那我只好用这个两个值进行一番太极拳，柔和在一起，得到一个估计值），要用这两个数据拟合出一个优化数据，就需要一个增益量Kg（权重系数），但是计算这个增益量很麻烦，用到上面提到的Q 、R，用来计算估计值误差（噪声）P和Kg。

其计算的基本原理大致如下图所示，左边的上一时刻的优化估计值分布模型（蓝色），和中间蓝色的预测值分布模型，以及右侧橙色的测量值分布模型，最终的目的是找到这一时刻的优化估计模型（灰色），这些分布都是正态分布，因此计算增益量其实就是找到一个优化正态分布模型的过程，这就要用到Q和R，计算出这个优化的正态模型的误差协方差P。

   基于以上两个基本公式，卡尔曼滤波最终要计算的是这一时刻的优化估计值。

X(k|k)=X(k|k-1)+Kg(k)*(Z(k)-Hx(k|k-1))            （3）

其中x(k|k)是这一时刻的优化估计值，x(k|k-1)是根据上一时刻的状态（优化后的值）对现在时刻状态的预测，在我们深度图例子中就是前面提到的观测者大脑根据上一时刻反馈出来的这一时刻的值（预测模型估计），Kg(k)就是我们要求的增益量（权重系数），Z（k）就是这一时刻设备测量值，即这一时刻的相机获取的深度图数据，H为测量系数，我们这个案例里面为1.

接下来那就算公式3中的每个值，以计算出x(k|k)，

根据1、2公式，可得

X(k|k-1)=AX(k-1|k-1)+Bu(k)+W(k)           （4）

本例 B=0  以及系统方程和观测方程计算不考虑噪声，所以 可变为

                              X(k|k-1)=AX(k-1|k-1)         

P(k|k-1)=A*P(k-1|k-1)A'+Q                      ( 5 )

P(k|k-1)为根据上一帧的优化模型估计的误差协方差P（k-1）来做预测的误差协方差，用来计算当前帧的增益量。

Kg(k)=P(k|k-1)*H'/(H*P(k|k-1)*H' + R).....  .( 6 )

这样通过公式4、5、6导入到3，就可以求出此时刻优化估计的值X(k|k)

那为了一个时刻P的预测，还需要计算当前时刻的P

p(k|k)=(I-Kg(k)*H)*P(k|k-1).....................( 5 )　　

总结起来就是下图的两个阶段（参考B站视频，注：表示符号与上面提到的公式不一样）

 

但在应用中，怎么去确定X、Z、A、H、Q、R呢？

X为我们想得到的状态结果，Z就是传感器的观测值

我自己最大的疑惑是观察方程里面怎么用状态方程里面的X表示，度娘很多文章，但是还是有疑问，但最终自己总结理解：Z(k)=HX(k)+V(k)    这个只是观察值的模拟表达式，观测方程是对传感器输出的描述，在直接观测模式下，H都是单位矩阵，z是通过传感器值输入的，只是在模拟里面或者没法直接测量（比如测火箭喷气温度，智能在接近的地方测）需要。

回到本深度图案例中，

X是我们深度图中的每个像素值和其与上一帧的变化值△x组成，A、B、U 、H取如下值

3.直方图均衡化

直方图均衡化(Histogram Equalization)是一种增强图像对比度(Image Contrast)的方法，其主要思想是将一副图像的直方图分布变成近似均匀分布，从而增强图像的对比度。直方图均衡化虽然只是数字图像处理(Digital Image Processing)里面的基本方法，但是其作用很强大，是一种很经典的算法。

其作用就是让灰蒙蒙的图，变得清晰些（灰度图），具体参考知乎直方图均衡化 - 知乎

使用的效果就像下图

 可使用Python库PIL中的ImageOps 的API实现

eq_img = ImageOps.equalize(img)