基于MATLAB的点云建筑物轮廓提取与基于平面探测法的点云建筑物提取

博客中轮廓提取使用的点云数据
建筑物平面检测使用的点云数据

**两个小的点云处理实验项目，（源码资源****有常（注意目前是有常哦）私我vx:xdsqczkyqs713
，第一个项目点云建筑物轮廓提取比较简单仅源码及相关测试文件20圆，第二个建筑物平面探测
源码30圆，带一份完整的方法实现论文40圆，报告仅供参考，可以就着写，节省时间，但宝儿你可千
万别傻乎乎的直接拿去上交鸭!）给的联系方式有常求源码的小可爱们加！非诚勿扰/<->认真脸）**

点云建筑物轮廓提取

1、txt文件点云数据导入：
使用uigetfile读取文件，加载后提取出点云x，y，z坐标，为了减少运算量提高速度，将坐标移动到原点，具体做法就是每个坐标值减去该坐标方向下的最小值：

[filename,pathname]=uigetfile({'*.txt';'*.ply';'*.las';},'载入点云数据');
if isequal(filename,0)||isequal(pathname,0)
    errordlg('没有选中文件','出错');
    return;
else
    file=[pathname,filename];
end
pointData = load(file);
x = pointData(:,1);
y = pointData(:,2);
z = pointData(:,3);
x = x-min(x);
y = y-min(y);
z = z-min(z);%移动到原点

2、点云图像散点图显示，主要是使用scatter3函数：

figure()
scatter3(x,y,z,'.')
grid on%显示网格
axis equal%每个坐标方向缩放比例设置为相同

3、依据高程着色，可以直接使用matlab自带的pcshow函数：

figure()
pcshow([x,y,z],'MarkerSize',500)%以点云文件格式显示，点大小设置为500，默认以z的大小作为点的颜色深度大小
axis off
saveas(gca,strcat(pathname,filename,'.png'));%保存影像化图片

4、将图片保存后，采用图像处理的方法来提取建筑物轮廓，主要是应用Canny算法的边缘检测：

BW = edge(I,'Canny');
figure()
imshow(BW)

5、在边缘检测的基础上，提取出最外围的轮廓即可

基于平面探测法的点云建筑物轮廓提取

1、方法流程图：

2、实验所采用的点云散点图显示结果如下：

3、使用pcshow函数，按照高程着色

4、基本上可以看出建筑物分布情况，大致可分为两个部分，按照高程切成两片区域：
在每个区域中，对区域进行切块遍历，采用平面探测的方法，检测平面点：

for  i=1:znum1%从z方向开始遍历,第一块，低洼区
    zind=find(z>zpiece1(i)&z<=zpiece1(i)+stepz);%先找到符合块内z方向的点
    if length(zind)>200%如果找到的点大于200，认为有可能拟合出平面
        for j=1:ynum1%再从y方向遍历
            yind=find(y>ypiece1(j)&y<=ypiece1(j)+stepy);%找到符合块内y方向的点
            yzCommon = intersect(yind,zind);%求交集，即得到既符合z方向又符合y方向的点
            if length(yind)>50&&length(yzCommon)>50%如果满足条件的点大于30个，认为有可能拟合出平面
                for k=1:xnum%最后从x方向开始遍历
                    xind=find(x>xpiece(k)&x<=xpiece(k)+stepx);%找到符合块内x方向的点
                    pointInd = intersect(yzCommon,xind);%求交集，找到符合块内xyz三个范围的点
                    if length(xind)>50&&length(pointInd)>50%如果符合的点数大于30，开始拟合
                        xt = x(pointInd);
                        yt = y(pointInd);
                        zt = z(pointInd);
                        [f,gof]= fit([xt,yt],zt,'poly11');%选择线性平面拟合
                        if gof.rmse<0.10%假如拟合的均方根误差小于0.15，认为拟合到了比较好的平面
                           planeNum = planeNum+1;%平面数加一
                           writematrix([xt yt zt],strcat('plane',num2str(planeNum)),'Delimiter',' ');%保存找到的平面点
                        else
                            continue%continue的意思是如果不符合条件，直接进入下一轮循环遍历
                        end
                    else
                        continue
                    end
                end
            else
                continue
            end
        end
    else
        continue
    end
end

5、最终检测的结果如下：