matlab对点云进行三角化(las格式读入)

1、LASread_and_convert.m文件（运行这个文件实现las读入为mat文件)

clear;
clc;
close all;
clear all;
A = LASreadAll('house_wall.las');%调用函数读取las文件
Point_City = [A.x,A.y,A.z];
save Point_City.mat;

 %将X1_buny赋值给一个p
    
    p = Point_City;
    
    %显示点云
    figure(1);
hold on
axis equal
title('点云','fontsize',14)
plot3(p(:,1),p(:,2),p(:,3),'g.')%这里颜色如果是 ‘r’就不行？？？
view(-37.5,30)


[t]=MyCrust(p);

%% plot of the oyput triangulation
figure(2)
        hold on
        title('三角化输出','fontsize',14)
        axis equal
        trimesh(t,p(:,1),p(:,2),p(:,3),'edgecolor','b')%plot della superficie trattata
        view(-37.5,30)

2、LASreadAll函数

function A = LASreadAll(infilename)
% LASREADALL reads in all variables from a LAS 1.x data file (used with lidar data) 
%
% INPUT
% infilename: input file name  (for example, 'myinfile.las') 
%           
% OUTPUT
% A:    This is a structure containing all the data in the file.
%         See the file documentation for more information:
% http://www.asprs.org/Committee-General/LASer-LAS-File-Format-Exchange-Activities.html
% 
% EXAMPLE
% A = LASreadAll('infile.las')
%
% Cici Alexander
%     September 2008 (updated 26.09.2008)
% Amy Farris (afarris@usgs.gov) 
%      November 2013 Substatially altered to read in all variables from the file

% This file has the capability to load data files using any 1 of 5 formats.
% However, I was only able to test it with one of the formats. 
% I included support for the rest of the file formats b/c I
% thought someone might find it useful.  I sure hope they work!!!
% If the code crashes, it may be because either the pointDataFormatID 
% value was not what I expected, or the value for pointDataRecordLength 
% is not what I think it should be based on the format ID.  
% If pointDataRecordLength is different than I expected, then you may need 
% to change the values in fseek, eg: (c+????)
% The number added to c should be the sum of bytes of all the variables 
% that occur before the variable currently being read in.
%
% Also, the file I tested this code with was LAS format 1.2, I 
% think this code will run on later versions of LAS.
%
% A brief explanation of the LAS file format:
% LAS files are binary, they begin with header information.  Included in
% the header is the size of the header (in bytes); this is called
% 'OffsetToPointData', refered to a 'c' in this code.  After c bytes the
% data begins with the first x value, then the first y value and so on...
% (exactly what data is included depends on the file format).  Then the
% file continues with the second x value, the second y value and so on...  
% The header tells you how many bytes of data there are for each data point
% ('pointDataRecordLength', also refered to as 'p' in this code).
% So to read in all the x values, you start at the beginig of the file and
% skip c bytes:  "(fseek(fid, c, 'bof');"  
% Then you read in one value and skip p-4 bytes 
% (each x value consists of 4 bytes) and then read the next x and so on:
% "X1 = fread(fid,inf,'int32',p-4);"
%
% This is my (Amy Farris') first attempt at reading in a binary file.  
% I depended strongly on Cici's orignal file (LASRead.m) at first.  Most 
% of the code after about line # 134 was written by me, as are these
% garrolous beginning comments.  I hope they help.

%% Open the file
fid =fopen(infilename);

% Check whether the file is valid
if fid == -1
    error('Error opening file')
end

%% Read in important information from the header
% Check whether the LAS format is 1.1
fseek(fid, 24, 'bof');
VersionMajor = fread(fid,1,'uchar');
VersionMinor = fread(fid,1,'uchar');
% afarris2011Aug20 changed the following line to read LAS1.2 files
% if VersionMajor ~= 1 || VersionMinor ~= 1 
if VersionMajor ~= 1  
    error('LAS format is not 1.*')
end

% Read in the offset to point data
fseek(fid, 96, 'bof');
OffsetToPointData = fread(fid,1,'uint32');

% Read in the point data fotmat ID
fseek(fid, 104, 'bof');
pointDataFormatID = fread(fid,1,'uchar');

% Read in the point data record length
fseek(fid, 105, 'bof');
pointDataRecordLength = fread(fid,1,'short');

% Read in the scale factors and offsets required to calculate the coordinates
fseek(fid, 131, 'bof');
XScaleFactor = fread(fid,1,'double');
YScaleFactor = fread(fid,1,'double');
ZScaleFactor = fread(fid,1,'double');
XOffset = fread(fid,1,'double');
YOffset = fread(fid,1,'double');
ZOffset = fread(fid,1,'double');

% The number of bytes from the beginning of the file to the first point record
% data field is used to access the attributes of the point data
c = OffsetToPointData;

% The number of bytes to skip after reading in each value is base