点云数据预处理&特征计算方法（全面总结整理）

1 点云数据分布特点

1.1 激光雷达点云特征

       （1）激光点云数据呈多层线性分布；
    （2）由于目标材质、颜色的不同，对激光点的反射强度也不同；
    （3）扫描距离越远，激光线束间隔越大，点云越稀疏；
    （4）用三维特征（x,y,z）表示点云数据；

1.2 4D毫米波雷达点云特征

        （1）用四维特征（R,θ,ψ,V）表示点云数据；（θ表示水平角，ψ表示俯仰角）

1.3 3D毫米波雷达点云特征

        （1）用三维特征（R,θ,V）表示点云数据；

1.4 基于深度图的三维点云

2点云数据预处理

2.1 筛选、过滤与分割

2.1.1 基于空间区域筛选

       建立感兴趣空间区域，将此空间以外的点剔除。图1 为三维激光雷达原始数据，坐标（0,0）为激光雷达所在位置，本文选取的高度范围为0.5~2 m，x 轴方向范围为0~25 m，y 轴方向范围为-10~10 m。筛选结果如图2所示，经筛选，剔除了树枝、地面等不在感兴趣区域内的障碍点，只保留了障碍物大致的轮廓，有利于提取车辆外围矩形。

2.1.2 基于网格法过滤与分割

2.2 点云数据下采样、上采样

2.2.1 点云数据下采样（稠密->稀疏）

       点云数据下采样的目的是为了将对全部点云的操作转换到下采样所得到的关键点上，从而达到降低计算量的目的。下采样的常用方法如下：

1）体素下采样

    体素下采样的原理如图1所示，首先将点云空间进行网格化，也称体素化，即图1(b)，网格化后的每一个格子称为体素，在这些划分为一个个极小的格子中包含一些点，然后对这些点取平均或加权平均得到一个点，以此来替代原来网格中所有的点，即图1©中蓝色的点。显然，网格选取越大则采样之后的点云越少，处理速度变快，但会对原先点云过度模糊，网格选取越小，则作用相反。

       体素下采样的特点是效率高，采样点分布相对比较均匀，同时可以通过控制网格尺寸控制点间距，但是不能精确控制采样点个数。
另外，对于二维数据下采样可采用网格下采样，原理类似。

2）均匀下采样

       均匀下采样有多种不同的采样方式，其中最远点采样是较为简单的一种，首先需要选取一个种子点，并设置一个内点集合，每次从点云中不属于内点的集合找出一点距离内点最远的点，如下图，这里的距离计算方式为该点至内点所有点的最小距离。这种方式的下采样点云分布均匀，但是算法复杂度较高效率低。

3）几何采样

       其原理是以点云的几何特征作为采样依据，这里以曲率为例。在点云中任意一点都存在某曲面，曲率计算示意图如图2所示，曲率越大，弧的弯曲程度越大，表示该地方的特征点越多，故在点云曲率越大的地方，采样点数越多，实现方法如下：首先计算每个点的K领域，然后计算点到领域点的法线夹角值，以此来近似达到曲率的效果并提高计算效率，因为曲率越大的地方，夹角值越大。设置一个角度阈值，当点的领域夹角值大于阈值时被认为是特征明显的区域，其余区域为不明显区域。对明显和不明显区域进行均匀采样，采样数分别为U*(1-V)和U*V,U是目标采样数，V是均匀采样性。其特点是计算效率高，且局部点云的采样是均匀的，同时稳定性高，使得采样结果的抗噪性更强。

4）随机下采样

       随机下采样的原理十分简单，如图所示，首先指定下采样的点数，然后进行随机点去除进行采样操作,得到图3(b)。随机下采样的特点是能控制输出点云的数量，但随机性太大，可能剔除点云的关键数据。

2.2.2 点云数据上采样（稀疏->稠密）

       点云数据上采样的目的是为了增加点云数量，以便更好的计算点云特征。

1）增采样

    增采样的原理如图4所示，当目前拥有的点云数据量较少时，如图4(a)，通过内插点云的方法对目前的点云数据对进行扩充，如图4(b)，达到保证基本形状不变的情况下增加点云。增采样的特点是可极大的增加点云数据，但由于内插点的不确定性会导致最后输出的结果不一定准确。

2）滑动最小二乘法采样

       滑动最小二乘法采样的原理是将点云进行了滑动最小二乘法的映射，使得输出的点云更加平滑。滑动最小二乘法的特点是适用于点云的光顺处理，但有时会牺牲表面拟合精度的代价来获得输出点云。
    具体的插值方法有很多，比如最近插值、线性插值等等，可以根据点云数据的特性及稀疏程度，以及算法实现的难易程度、耗时等，来设计实用性的增采样方法。

2.2.3 总结

       在下采样方法中，以体素化网格采样方法最为常用，因为其速度快，代码量少，且满足大多数时的点云处理要求；均匀采样虽然精度高，当耗时高，可以用于更追求精度的场合下；几何采样由于使用不多，方法很多，这里只是简答介绍了一下曲率采样，比较适用于不规则的且丰富表面特征的点云数据计算；随机下采样由于能准确控制点云的输出数量，但过于随机，较少使用；
    增采样用于增加点云数据，更适合用于解决曲面重建时点云数量缺少的问题；而滑动最小二乘法同样是对点云数量的扩充，但主要是对点云形状进行平滑处理，所以更适合用来对点云结构进行优化。

2.3 点云滤波（去除离群点）

2.3.1 半径滤波器

       设定某个点半径范围，其邻点数量少于设定的阈值，就认为该点为离群点。

2.3.2 统计滤波器

       稀疏离群点去除，基于计算输入数据中的点与邻域距离的分布。对于每个点，计算它与所有邻居之间的平均距离。假设结果分布是具有平均数和标准差的高斯分布，所有平均距离在由全局距离平均数和标准差定义的区间之外的点都可以被视为离群点，并从数据集中修剪掉。

2.3.3 条件滤波器

       可以设定删掉符合指定条件的点。

2.3.4 直通滤波器

       直通滤波器，就是在点云的指定维度上设置一个阈值范围，将这个维度上的数据分为在阈基于高斯核的卷积滤波实现 高斯滤波相当于一个具有平滑性能的低通滤波器。

2.3.5 双边滤波

       双边滤波算法，是通过取邻近采样点的加权平均来修正当前采样点的位置，从而达到滤波效果。同时也会有选择地剔除部分于当前采样点“差异”较大的相邻采样点，从而达到保持原特征的目的。

2.3.6 移动最小二乘法光滑滤波

       根据点云，找到一个函数，这个函数由一系列的局部的函数组成（这个工作是通过移动最小二乘法完成的），这个函数的特性是光滑的，然后将点投影到函数上，就完成了平滑的目的。

2.4 点云特征提取

2.4.1 点云聚类方法

       1）基于划分的聚类方法就是将一个包含有 N 个数据对象的数据集划分为 K 个组，每一组代表一个聚类且 K≤N.
    2）基于层次的聚类方法就是将待处理的数据集不断分解并创建层次；
    3）基于密度的方法，根据数据对象的密集程度进行聚类，当密度超过一定阈值时会停止，这种方法可以发现不规则形状的数据聚类组；
    4）基于网格的聚类方法是将数据集所在的空间划分为由单元格组成的网络结构，由单元格取代单个数据对象。此种方法的优势在于面对数据量较大的数据集时可以进行快速聚类节省时间。

DBSCAN算法

       在x轴方向依据点与点的距离差值聚类，在此基础上再在y 轴方向根据点与点的距离差值聚类。最后输出障碍点数聚类集合m={m1,m2,m3,…,mn}。
    然后对输出数据集中的每个障碍点集合分别应用DBSCAN聚类算法，对障碍点进行进一步分类.DBSCAN聚类算法需要事先输入过滤噪声的核心点领域距离阈值Eps 以及核心点邻域范围内最少点个数的阈值MinPts:

2.4.2 点云包围框提取

       基于点云聚类的结果，对聚类点进行包围框提取，常见的包围框有圆形包围框、矩形包围框、最小凸包围框，提取方法如下：

1）圆形包围框

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function: cal Circular bounding box
% Author: ZZ
% Create Time: 2022.9.21
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [cenp,r] = calCircularBoundingbox()
%     Ritter's算法如下：
%     1.从点集中随机选出两个点作为直径对圆进行初始化。
%     2.判断下一个点p是否在圆中，如果在则继续本步骤，如果不在则进行步骤3。
%     3.使用p作为新圆的一个边界点，另一个边界点为距离p最远的圆上的点，使用这两个点作为直径构造新圆。
%     4.继续步骤2，直到遍历完所有点。
%     返回值cenp中心点位置,r半径

    n = 100;
    p = rand(n,2);

    p1 = p(1,:);
    p2 = p(2,:);
    r = sqrt((p1(1)-p2(1))^2+(p1(2)-p2(2))^2)/2;
    cenp = (p1+p2)/2;

    for i = 3:n
        newp = p(i,:);    
        d = sqrt((cenp(1)-newp(1))^2+(cenp(2)-newp(2))^2);  
        if d > r
            r = (r+d)/2; % 该算法求出的圆比最优圆大概会大个5%到20%左右
            cenp = cenp + (d-r)/d*(newp-cenp);
        end    
    end
    figure(9)
    clf
    hold on;
    plot(p(:,1),p(:,2),'o');
    x0=cenp(1);
    y0=cenp(2);
    theta=0:0.01:2*pi;
    x=x0+r*cos(theta);
    y=y0+r*sin(theta);
    plot(x,y,'-',x0,y0,'.');
    axis equal
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

2）矩形包围框

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function: cal Rectangular bounding box
% Author: ZZ
% Create Time: 2022.9.22
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Rectangular] = calRectangularBoundingbox()
    n = 50;
%     rng('default');
    p = rand(n,2);

    ind = convhull(p(:,1),p(:,2)); % 系统函数求凸包点
    hullNum = length(ind);
    hull = p(ind,:);               % 随机点凸包
    area = inf;
    for i = 2:hullNum
        p1 = hull(i-1,:);     %凸包上两个点
        p2 = hull(i,:);
        k1 = (p1(2)-p2(2))/(p1(1)-p2(1));     %连接两点的直线，作为矩形的一条边
        b1 = p1(2)-k1*p1(1);

        d = abs(hull(:,1)*k1-hull(:,2)+b1)/sqrt(k1^2+1);  %所有凸包上的点到k1,b1直线的距离
        [h ind]=max(d);                     %得到距离最大的点距离，即为高，同时得到该点坐标
        b2 = hull(ind,2)-k1*hull(ind,1);    %相对k1,b1直线相对的另一条平行边k1,b2;
        k2 = -1/k1;                         %以求得的直线的垂线斜率

        b = hull(:,2)-k2*hull(:,1);         %过凸包所有点构成的k2,b直线系
        x1 = -(b1-b)/(k1-k2);               %凸包上所有点在已求得的第一条边的投影
        y1 = -(-b*k1+b1*k2)/(k1-k2);

        x2 = -(b2-b)/(k1-k2);               %凸包上所有点在已求得的第二条边的投影
        y2 = -(-b*k1+b2*k2)/(k1-k2);

        [junk indmax1]=max(x1);             %投影在第一条边上x方向最大与最小值
        [junk indmin1]=min(x1);

        [junk indmax2]=max(x2);             %投影在第二条边上x方向最大与最小值
        [junk indmin2]=min(x2);

        w = sqrt((x1(indmax1)-x1(indmin1))^2+(y1(indmax1)-y1(indmin1))^2);    %矩形的宽
        if area >= h*w                      %使面积最小
            area = h*w;
            pbar=[x1(indmax1) y1(indmax1);  %矩形四个角点
                  x2(indmax2) y2(indmax2);        
                  x2(indmin2) y2(indmin2);
                  x1(indmin1) y1(indmin1)];            
        end
    end
    Rectangular = pbar; % 返回矩形框顶点坐标
    pbar(5,:) = pbar(1,:);
    figure(4)
    clf
    plot(p(:,1),p(:,2),'.');hold on;
    plot(hull(:,1),hull(:,2),'bo');hold on;
    plot(pbar(:,1),pbar(:,2),'r')
    axis equal;
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

3）最小凸包围框

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function: cal Minimum convex bounding box
% Author: ZZ
% Create Time: 2022.9.22
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hull] = calMinConvexBoundingbox()
    % Graham扫描法计算最小凸包围框
    % 找到离散点中，保证y坐标最大的情况下，x坐标最小的点，记做A点
    % 以A点为原点，x轴正反向射线顺时针扫描，找到旋转角最小时扫描到的点，记做B点，循环
    n_points = 25; % total number of points
    rng('default');
    points_source = rand(n_points,2);
    points = points_source;
    
    figure(1);clf;
    plot(points(:,1),points(:,2),'.k');
    hold on;
        
    hull = [];
    inner = [];
    undefined = (1 : n_points)';
    % As the origin point of the convex hull we take point with the lowest Y value.
    % Origin point is a priori lies on the border of the convex hull.
    [~, lowest_point_id] = min(points(:, 2));

    % Calculate angles between X axis and each segment, connecting the origin
    % with all of the points.
    u = points(lowest_point_id,1:2);
    angles = zeros(n_points,1);
    for i = 1 : n_points
        v = points(i, 1 : 2);    
        w = v - u;
        angles(i) = atan2d(w(2), w(1));
    end
    
    % Sort points in according to it's angles with the X axis.
    undefined = [undefined,angles,points(:,1)]; % [序号，夹角，x坐标值]
    undefined = sortrows(undefined, 2);         % 按照角度排序
    undefined = undefined(:, 1);                % 取排序后的序号

    % First point will be on the convex hull.
    hull = undefined(1);                        % 凸包点
    plot(points(hull(end),1),points(hull(end),2),'r*');
    text(points(hull(end),1),points(hull(end),2),num2str(hull(end)));
    undefined = [undefined; undefined(1)];
    undefined(1) = [];
    
    % Loop through all points and remove those of them which are not on the border
    % of the convex hull.
    h1 = [];count = 0;
    while length(undefined) >= 2
        while true
            count = count + 1; % 循环次数
            test_ids = [hull(end);undefined(1:2)];
            test_points = points(test_ids,:);
            h1 = plot(test_points(:,1),test_points(:,2),'bo');
            if isCounterclockwise(test_points(1,:),test_points(2,:),test_points(3,:))
                break
            end
            inner = cat(1,inner,undefined(1)); % 内点
            if ~isempty(h2)
                delete([h2(hull(end)),h3(hull(end))])
            end
            undefined(1) = [];
            undefined = cat(1, hull(end), undefined);
            hull(end,:) = [];
            if ~isempty(h1)
                delete(h1)
            end
        end
        if ~isempty(h1)
            delete(h1)
        end
        hull = cat(1,hull,undefined(1));
        undefined(1) = [];
        h2(hull(end)) = plot(points(hull(end),1),points(hull(end),2),'r*');
        h3(hull(end)) = text(points(hull(end),1),points(hull(end),2),num2str(hull(end)));
        drawnow
    end
    hull = cat(1,hull,undefined(1));
    plot(points(hull,1),points(hull,2),'.-r');
    count

    function result = isCounterclockwise(a, b, c) % 判断是否逆时针方向
        result = ((b(1) - a(1)) * (c(2) - a(2)) - (b(2) - a(2)) * (c(1) - a(1)) >= 0);
    end

end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

2.4.3 点云直线特征提取

       直线特征提取方法有：最小二乘拟合直线提取、Hough变换直线提取等。Hough变换直线提取的代码实现方法如下：

y = 1:0.5:100;
y = y + randn();
x = 3 .* y + 14 + randn();
figure(11)
clf
plot(x,y,'*')
findLineByHoughTransform(x,y);

function [R,T] = findLineByHoughTransform(x,y)

% R = x * cosT + y * sinT

        Rmax = round(sqrt(max(x)^2 + max(y)^2));
        
        % Hough Transform step
        thetaStep = 1; % 角度°
        rangeStep = 1;
        theta = -90 : thetaStep : 90;
        range = -Rmax : rangeStep : Rmax;
        
        % 初始化参数空间矩阵
        A = zeros(size(theta,2),size(range,2));
        xMax = zeros(size(theta,2),size(range,2));
        xMin = zeros(size(theta,2),size(range,2));

        for it = 1:size(theta,2)
            for ip = 1:length(x)
                r = x(ip) *cosd(theta(it)) + y(ip) *sind(theta(it));
                Rindex =  round((r + Rmax)/rangeStep) + 1;
                A(it,Rindex) = A(it,Rindex)+1;
                if A(it,Rindex) == 1
                    xMax(it,Rindex) =  x(ip);
                    xMin(it,Rindex) =  x(ip);
                else
                    xMax(it,Rindex) =  max(x(ip),xMax(it,Rindex));
                    xMin(it,Rindex) =  min(x(ip),xMin(it,Rindex));
                end
            end
        end
        
        % 找到通过直线最多的参数点
        [indexT,indexR] = find(A == max(max(A)));
        R = range(indexR);
        T = theta(indexT);
        
         figure(11)
         hold on
         for it = 1:length(indexT)        
            k = -1/tand(theta(indexT(it)));
            b = range(indexR(it))/sind(theta(indexT(it)));
            x = xMin(indexT(it),indexR(it)):0.1:xMax(indexT(it),indexR(it));
            y = k.*x + b;
            plot(x,y,'ro')
            hold on;
         end


end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

2.4.4 聚类点云的分布特征

       1）点云离散距离主要表征边缘点云偏离程度；
    2）三维协方差主要表征点云在各方向上的相关性，而特征值表征该方向的权重；
    3）惯性张量矩阵主要表征点云整体分布的稳定性，可用于消除噪声影响；

       其他相关特征计算：
https://blog.csdn.net/jane0819/article/details/130983616