RGBD相机实用问题

奥比中光3D视觉开发者社区

于 2022-09-16 11:00:00 发布

阅读量1.6k

点赞数 3

文章标签：图像处理计算机视觉 opencv

原文链接：https://developer.orbbec.com.cn/v/blog_detail?id=777

版权

作者：玉林
编辑：3D视觉开发者社区

文章目录

1.深度图像格式和存储
2.查看16bit png格式深度图像的方法
3.相机模型
4.相机标定参数
5.内参的几何意义
6.图像裁剪或旋转后，内参如何相应改变
7.结构光深度相机的测量精度
8.深度相机的深度图像和RGB图像对齐
9.深度点云平面拟合

经常看到论坛上有初学者针对奥比中光的RGBD相机提出问题，这篇博客试图回答一些常见问题。

原文地址在奥比中光的3D视觉开发者社区：https://developer.orbbec.com.cn/v/blog_detail?id=777，转载请注明出处。

1.深度图像格式和存储

深度图像是一个单通道16bit二维无符整形数组，保存的是三维空间点的深度 $Z$ 值，深度单位通常是毫米。在C++中用 unsigned short*表示数组对应内存，内存长度为图像尺寸W×H×2个Byte，所以640*480的深度图的长度为600KB。

深度图可用C++函数fwrite和fread保存和读取raw data，但比较方便的方法是用单通道png格式存储图片(是一张灰度图)。因为png可以存储每像素16bit（或2个Byte）的图像，而且png图像是无损压缩的，这点对深度图来说非常重要。另外，png作为通用格式可以用图像浏览软件进行查看。

Matlab的imread和imwrite函数对可以对png文件进行读写。OpenCV可以把16bit无符整型数组放在Mat结构里，调用imwrite函数进行png图像格式存储:

cv::Mat cv_depth;
uint16_t* pDepth;
//......
cv_depth = cv::Mat(hei, wid, CV_16UC1, pDepth);
cv::imwrite("4.png", cv_depth);

需要注意的是用OpenCV读取16bit png深度图时，要使用参数CV_LOAD_IMAGE_UNCHANGED 或 -1：

CV::Mat cv_depth = imread("4.png",CV_LOAD_IMAGE_UNCHANGED);

这样图像是CV_16UC1类型数据，可以用cv_depth.at(i,j)访问像素。

2.查看16bit png格式深度图像的方法

尽管png格式深度图可以用“画图”等工具打开，但是可能是下面的乌漆麻黑的效果：

在这里插入图片描述
原因是16bit能表达的范围是 $0,2^{16}-1]$ ，深度图最大值通常只有10000(毫米)，所以需要把深度范围拉伸到合适的范围内显示。比如把一幅深度图的最小和最大值之间的数值线性通过变换拉伸到[0,255]之间。

在Matlab中就比较简单，imshow中使用"[]"就能做显示范围的拉伸：

>> dp=imread('4.png');
>> figure,imshow(dp,[]);

结果如下边左图所示，而一些图像浏览软件如ImageJ.exe可以做类似的操作，如下边右图。

在这里插入图片描述
而下面的代码可以用一个映射表map，把0~3000范围内的深度值变换为伪彩色的形式：

>> map=jet(3000);
>> figure,imshow(dp,map);

这样，图中冷暖色调的变化对应深度由小到大，距离由近及远的变化（蓝色背景无深度，深度值为0）：
在这里插入图片描述

3.相机模型

相机绕着物体拍摄，相机坐标系和物体坐标系有相对的位姿变化，6个自由度的RT，包括3个旋转和3个平移。成像过程光线通过透镜中心投影在成像面(想像成相机底片或CCD)上成像，一般为了方便起见把成像面移到物方同侧，如下边右图所示：
在这里插入图片描述
图中相机坐标系xyz轴用RGB颜色表示，深度相机中的深度值就是空间点的zz坐标。成像面是一个2D图像，坐标系原点在图像左上角。成像面到相机坐标系中心的距离称为焦距。

4.相机标定参数

相机内参如果不考虑畸变通常是4个数，即焦距 $f_x,f_y$ 和光心坐标 $(u, v)$ （内参也可以是5个数，一般3个就够了，即 $f_x=f_y$ ）组成的3*3的内参矩阵：

$K=\left[\begin{matrix} f_x & 0 & u \\ 0 & f_y & v\\ 0 & 0 & 1\end{matrix}\right]$

内参矩阵提供了“物点-像点-光心o”三者共线关系的一种表达，内参要通过标定过程进行估计。经过标定的相机每个像素对应了一条相机坐标系原点发出的射线，相机就变成了一个三维光学测量工具。

内参可用于3D空间点 $(X, Y, Z)$ 到2D图像点 $(x, y)$ 的投影：
$\lambda[x,y,1]^T=K[X,Y,Z]^T$
具体来说：
$x=f_x*X/Z+u\\ y=f_y*Y/Z+v$
可见相机成像是个降维过程，丢失了三维空间信息，但是深度相机把深度 $Z$ 保留在了 $(x, y)$ 坐标位置。
对于深度相机，已知深度图上像素位置 $(x, y)$ 的深度 $Z$ 值，亦可利用内参反求出空间点的3D坐标 $(X, Y, Z)$ ：
$X=(x-u)* Z/f_x\\ Y=(y-v)* Z/f_y\\ Z=Z$
奥比中光生产的每台RGBD相机，出厂时都做过标定，标定参数都存储在设备里。对于开发者来说，内参可以调用SDK直接从设备里读取，不用自己重新标定，除非你有更高的精度需求。

5.内参的几何意义

通过相机模型四棱锥的俯视和侧视图，可以看出内参和图像分辨率之间的几何关系：
在这里插入图片描述

W和H分别为图像的宽和高，u和v是图像中光心在图像上的坐标，焦距fx和fy表示原点o到相平面的距离。内参中的数字的单位是像素，u和v通常大约是图像分辨率W和H的一半。

通过W与fx的关系以及H与fy的关系，可以估计出相机的水平视角HFOV和垂直视角VFOV。

例子某台AstraPlus RGBD设备的深度图像的内参读出来数值如下：

$f_x=f_y=540.393,u=320.617,v=228.87$

深度图像分辨率为W=640,H=480（一般也称为VGA格式）。可以看出u和v大约为W和H的一半。

如果不知道图像分辨率是VGA还是QVGA（W:H=320:240），或者HD和FullHD等常见格式中的哪一种，根据u和v也可以推测出来。

此时FOV计算结果为：

HFOV=2*atan(320/540.393)180/pi=61.2648
VFOV=2atan(240/540.393)*180/pi=47.8940

6.图像裁剪或旋转后，内参如何相应改变

相机一般在某个分辨率下进行内参标定（标定过程需要对着标定板采集照片），如果在使用时需要对图像进行诸如缩放、旋转和裁剪等操作，内参矩阵中的4个数也需要做相应的修改。这样修改后的内参作用在深度图像上才能恢复正确的空间点坐标。

实际上，根据前边的内参的几何意义一节，可以看出内参矩阵中的焦距 $f_x,f_y$ 和光心坐标 $(u, v)$ 应该如何变化。

缩放：图像尺寸W和H发生变化，相机模型空间四棱锥发生整体缩放，焦距 $f_x, f_y$ 和光心坐标 $(u, v)$ 要经过相应的等比例变换，此时视角不变。假设图像缩放比例为 $s$ :

$K=\left[\begin{matrix} sf_x & 0 & su \\ 0 & sf_y & sv\\ 0 & 0 & 1\end{matrix}\right]$

±90°或180°旋转：右转90°，焦距的两个值要交换位置，

$K=\left[\begin{matrix} f_y & 0 & H-v \\ 0 & f_x & u\\ 0 & 0 & 1\end{matrix}\right]$

左转90°，

$K=\left[\begin{matrix} f_y & 0 & v \\ 0 & f_x & W-u\\ 0 & 0 & 1\end{matrix}\right]$

右转180°，

$K=\left[\begin{matrix} f_x & 0 & W-u \\ 0 & f_y & H-v\\ 0 & 0 & 1\end{matrix}\right]$

垂直翻转： $f_x,f_y$ 不变，新的光心坐标 $(u ’, v^{'}) = (W - u, v)$
$K=\left[\begin{matrix} f_x & 0 & W-u \\ 0 & f_y & v\\ 0 & 0 & 1\end{matrix}\right]$

裁剪：新图像尺寸为 $(W^{'}, H^{'})$ ,焦距 $f_x,f_y$ 不变，光心坐标 $(u, v)$ 在新图像中的位置为 $(u - d x, v - d y)$ :
$K=\left[\begin{matrix} f_x & 0 & u-dx \\ 0 & f_y & v-dy\\ 0 & 0 & 1\end{matrix}\right]$

这时由于裁剪，相机的视野(FOV)变窄了。

总结在下图中：
在这里插入图片描述

7.结构光深度相机的测量精度

结构光深度相机利用三角测距原理，可以看成是个双目相机，只是左眼换成了激光器发射的激光散斑。双目估计深度的原理是：估计左眼的像素点在右眼中的偏移，称为视差(disparity)，逆深度值的变化量 $\Delta1/Z$ 和视差 $d X$ 成正比。如下图， $Z_2$ 是已知参考图(正对着一个平面采集的)上的已知深度，通过散斑视差求深度 $Z_1$ :
在这里插入图片描述

假设视差 $d X$ 只能估计到1/8像素精度，现在要反推在不同 $Z_1$ 位置它对应的变化量，即感知精度。公式中左右眼距离即基线B=75mm，图像焦距f=555pixel@VGA分辨率，已知 $Z_2$ 求 $Z_1$ ，公式：

$Z_2=uint16(1/(1/Z_1-0.125/(75*555)))$ 。根据 $Z_1$ 位置的不同，1/8像素视差对应的深度值的增量如下：

$Z_1=2000$ 时， $Z_2=2012$ ，增量是 $12 mm$ ;
$Z_1=4000$ 时， $Z_2=4049$ ，增量是 $49 mm$ ;
$Z_1=8000$ 时， $Z_2=8197$ ，增量是 $197 mm$ 。

可见，双目相机测距的误差，随着距离增加不是线性增加的，距离越远增大得越快，即相机的感知精度越差。实际上，结构光深度相机虽然有0到10000mm的量程，能输出的有效整数深度值也只不过有800多个。不信你可以把结构光深度相机所能输出的深度值统计一下试试看。

这样看结构光深度相机得到的整数深度值分辨率只有在近距时是 $1 mm$ 左右精度的，如果拍摄距离较远的光滑曲面得到的可能是阶梯状表面。这就需要对深度图小心翼翼地做一些滤波处理，比如双边滤波(Bilateral filtering)，以保证既能得到光滑的表面又不至于丢失过多的深度细节。

8.深度相机的深度图像和RGB图像对齐

深度相机的深度图像来自红外摄像头，另外设备还有RGB传感器用于点云着色。比如奥比中光的Astra相机结构如下：

在这里插入图片描述

可以看出红外和RGB摄像头二者有安装距离，所以从实际采集的图像上看存在视差。为了把深度和RGB图像对齐，消除视差，得到一个真正的RGB-D图像，即物体表面的颜色和它的深度在二维图像上像素级精确对应，需要做一个转换：首先把红外相机坐标系下的深度先换为空间点云，再通过刚性变换转换到RGB摄像头的坐标系，并最终投影到RGB图像的二维图像坐标系，形成一张在RGB相机坐标系下的深度图。

RGB和深度摄像头坐标系之间的刚性变换要用到二者之间的外参即RT关系，这也是经过出厂标定保存在设备里的。开发者可以调用SDK读取这个外参。另外，调用SDK读取RGBD序列时，通过函数可以确定是否开启“深度-彩色图对齐(D2C)”功能，它是用设备的芯片计算前面描述的过程。

如果你用的是OpenNI的SDK，方法如下：

openni::Device*       pDevice;
pDevice = new openni::Device();
...
pDevice>setImageRegistrationMode(openni::IMAGE_REGISTRATION_DEPTH_TO_COLOR);

openni::VideoFrameRef frame;
int rc = pDepthStream->readFrame(&frame);
//这时得到的深度图就是和RGB对齐的
uint16_t * pDepth = (uint16_t*)frame.getData();

9.深度点云平面拟合

RANSAC是机器视觉领域中为大家熟知的算法，以至于点云拟合平面首先被想到的就是RANSAC。实际上没有任何先验的RANSAC可能拟合出不正确的结果。当平面空间位置大致已知的时候，也可以利用鲁棒最小二乘法拟合平面。它是一种带权重的最小二乘迭代方法，每个点的权重和它的拟合误差有关，距离远的权重小。它的想法是通过改变权重抑制那些离群点(outlier)对拟合平面的影响。比如，权重和距离x的关系可以用 $w=1/(x^2+600)$ ，把这个关系画出来看看：

x = [-300:300];
w = 1./(x.^2+600);
w = w/sum(w);
figure,plot(x,w);
xlabel('distance'),ylabel('weight');

长这样：

在这里插入图片描述

在超过一定距离时，权重衰减得很快，但也不完全为0。

深度点云 $(X, Y, Z)$ 满足的平面条件可以写成 $Z = a X + bY + c$ ，但是对于平面和光轴接近垂直的时候导致无法准确拟合平面，也就是说这个方程表达对平面法线朝向要求苛刻，所以一般把平面的表示为 $A X + B Y + CZ + D = 0$ 。平面有3个自由度，平面法线为单位法线矢量： $\vec n=(A,B,C)^T$ ，包含有2个自由度，另1个自由度是光心到平面距离DD。光心到平面上最近点的坐标是 $(- D A, - D B, - D C)$ 。

平面拟合过程需要求解方程 $X,Y,Z,1][A,B,C,D]^T=0$ 。把点云凑成一个n×4的矩阵 $P$ ，每行一个齐次坐标表示的三维点，行数n是参与平面拟合的点的个数。这样方程组变成 $PX = 0$ 。

求解形如 $A X = 0$ 的齐次线性方程组要用SVD求解，矩阵A的最小Eigenvalue对应的本征向量Eigenvector，又称为A的null space。它的几何解释：找到一个向量尽可能和A中的每一列都正交(垂直)。这是《计算机视觉中的多视觉几何》一书最重要的知识点之一。实际上是对 $A^TAX=0$ 求解，因为 $A^TA$ 是4维矩阵，SVD求解更容易。

用前面的人体模型深度图（仿真合成的深度）进行地面拟合实验：

dp = imread('data\4.png');dp = single(dp);
Kc=[520,0,240;
    0,520,320;
    0,0,1];

[H,W,~]=size(dp);

% depth to 3D points
[x2d,y2d] = meshgrid(1:W,1:H);
X = (x2d-Kc(1,3)).*dp/Kc(1,1);
Y = (y2d-Kc(2,3)).*dp/Kc(2,2);
 
% 已知平面大致位置
ABCD = [0,-1,0,800];

% 通过距离选择包含平面的点进行拟合
pt = [X(:),Y(:),dp(:)];
dis = ABCD(1)*pt(:,1)+ABCD(2)*pt(:,2)+ABCD(3)*pt(:,3)+ABCD(4);
id = abs(dis)<400&dp(:)>0;
pt = pt(id,:);

% weighted MSE
wei = ones(size(pt,1),1);
for round=1:20
    ABCD = fit_plane(pt,wei);
    dis = ABCD(1)*pt(:,1)+ABCD(2)*pt(:,2)+ABCD(3)*pt(:,3)+ABCD(4);
    wei = 1./(dis.^2+600); wei = 1*wei/sum(wei);
end

% 标记平面点，选距离1CM内的点
id=abs(ABCD(1)*pt(:,1)+ABCD(2)*pt(:,2)+ABCD(3)*pt(:,3)+ABCD(4))<10;
figure,hold on
plot3(pt(:,1),pt(:,2),pt(:,3),'.'),axis equal
plot3(pt(id,1),pt(id,2),pt(id,3),'g.')
view(0,-70)

%% fit plane with SVD
function ABCD = fit_plane(pt,wei)
pnum=size(pt,1);
A=[pt,ones(pnum,1)];
W=repmat(wei,1,4);
[u,~,~]=svd(A'*(W.*A));
%last column of u
ABCD = u(:,4);
ABCD = ABCD/norm(ABCD(1:3));
end