function lasHeaderInfo = read_las_header(lasFilename)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Reads and displays the header information from a lidar data file in    %
%  the American Society for Photogrammetry & Remote Sensing (ASPRS) LAS   %
%  format, as specified in the following document: "LAS Specification     %
%  Version 1.2" ASPRS, April 29, 2008.  It is backwards-compatible,       %
%  with respect to the LAS version (i.e., header information from lidar   %
%  files in LAS 1.0 and 1.1 formats will display fine).  Requires the     %
%  file "coordinate_sys.mat" to display coordinate system info.           %
%                                                                         %
%  Input:                                                                 %
%     lasFilename: Path and filename of LAS file                          %
%                  Example of calling this function:                      %
%                  headerInfo = read_las_header('C:\my_lidar_file.las');  %
%                                                                         %
%  Output:                                                                %
%     lasHeaderInfo: Up to 40x2 cell arrray containing all the info read  %
%                    in from the LAS file header. If coordinate system    %
%                    info is missing, then lasHeaderInfo will be a 36x2   %
%                    cell array.                                          %
%                                                                         %
%  Chris Parrish                                                          %
%  Created:  5/13/2008                                                    %
%  Modified: 2/26/2009: Minor edits in preparation to post on MATLAB      %
%                       Central. Tested in MATLAB v7.0.1.15 (R14) and     %
%                       Octave v3.0.3.                                    %
%            9/2/2011:  Fixed bug that caused Version Major, Version      %
%                       Minor, and Point Data Format IDs to not be        %
%                       read in correctly. Also now displays vertical     %
%                       datum info and units, if polulated in             %
%                       GeoKeyDirectoryTag.                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Open the LAS file
fid=fopen(lasFilename);
if fid == -1 
    error('LAS file could not be opened; check filename and path.')
end

% Read the header info. Numeric data is always read in as a double because 
% built-in MATLAB funtions typically expect doubles.
lasHeaderInfo(1,1:2) = [{'File Signature'}; ...
    {sscanf(char(fread(fid,4,'char*1=>char*1')'),'%c')}];
lasHeaderInfo(2,1:2) = [{'File Source ID'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(3,1:2) = [{'Global Encoding'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(4,1:2) = [{'Project ID - GUID data 1'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(5,1:2) = [{'Project ID - GUID data 2'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(6,1:2) = [{'Project ID - GUID data 3'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(7,1:2) = [{'Project ID - GUID data 4'}; ...
    {sscanf(char(fread(fid,8,'char*1=>char*1')'),'%c')}];
lasHeaderInfo(8,1:2) = [{'Version Major'}; ...
    {fread(fid,1,'uint8=>double')}];
lasHeaderInfo(9,1:2) = [{'Version Minor'}; ...
    {fread(fid,1,'uint8=>double')}];
lasHeaderInfo(10,1:2) = [{'System Identifier'}; ...
    {sscanf(char(fread(fid,32,'char*1=>char*1')'),'%c')}];
lasHeaderInfo(11,1:2) = [{'Generating Software'}; ...
    {sscanf(char(fread(fid,32,'char*1=>char*1')'),'%c')}];
lasHeaderInfo(12,1:2) = [{'File Creation Day of Year'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(13,1:2) = [{'File Creation Year'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(14,1:2) = [{'Header Size'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(15,1:2) = [{'Offset to point data'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(16,1:2) = [{'Number of variable length records'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(17,1:2) = [{'Point Data Format ID'}; ...
    {fread(fid,1,'uint8=>double')}];
lasHeaderInfo(18,1:2) = [{'Point Data Record Length'}; ...
    {fread(fid,1,'uint16=>double')}];
lasHeaderInfo(19,1:2) = [{'Number of point records'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(20,1:2) = [{'Number of points return 1'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(21,1:2) = [{'Number of points return 2'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(22,1:2) = [{'Number of points return 3'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(23,1:2) = [{'Number of points return 4'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(24,1:2) = [{'Number of points return 5'}; ...
    {fread(fid,1,'uint32=>double')}];
lasHeaderInfo(25,1:2) = [{'X scale factor'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(26,1:2) = [{'Y scale factor'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(27,1:2) = [{'Z scale factor'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(28,1:2) = [{'X offset'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(29,1:2) = [{'Y offset'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(30,1:2) = [{'Z offset'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(31,1:2) = [{'Max X'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(32,1:2) = [{'Min X'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(33,1:2) = [{'Max Y'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(34,1:2) = [{'Min Y'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(35,1:2) = [{'Max Z'}; ...
    {fread(fid,1,'double=>double')}];
lasHeaderInfo(36,1:2) = [{'Min Z'}; ...
    {fread(fid,1,'double=>double')}];

% Read the coordinate system info from the variable length records, if 
% available
filePos = lasHeaderInfo{14,2};  % Set file position based on header size
nVarRec = lasHeaderInfo{16,2};  % Number of variable length records
if ~(nVarRec > 0) || ~exist('coordinate_sys.mat')
    disp('Note: Cannot display coordinate system info because')
    disp('coordinate_sys.mat file was not found, or the LAS file does')
    disp('not contain any varaible length records.')
else
    load('-mat', 'coordinate_sys')
    for i = 1:nVarRec
        fseek(fid, filePos, 'bof');
        rsvd = fread(fid,1,'uint16=>double');
        userId = sscanf(char(fread(fid,16,'char*1=>char*1')'),'%c');
        recId = fread(fid,1,'uint16=>double');
        lenRec = fread(fid,1,'uint16=>double');
        desc = sscanf(char(fread(fid,32,'char*1=>char*1')'),'%c');
        filePos = filePos + lenRec + 54;
        % If the current variable length record contains the
        % GeoKeyDirectoryTag, then display the georeferencing info. Note:
        % GeoKeyDirectoryTag = 34735.
        if strcmp(userId,'LASF_Projection') && recId == 34735
            nTags = floor(lenRec/2);
            % Jump to the correct position in the file
            newFilePos = filePos - lenRec;
            fseek(fid, newFilePos, 'bof');
            nCodesFound = 0;     % Number of codes found that match tags
            for j = 1:nTags
                tag = fread(fid,1,'uint16=>double');
                % Read through the CS and units codes and see if the tag
                % matches a known code; if so store it
                index = find(projectedCsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Projected Coordinate System'}; ...
                        projectedCsNames(index)];
                end
                index = find(geographicCsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Geographic Coordinate System'}; ...
                        geographicCsNames(index)];
                end
                index = find(vertCsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Vertical Coordinate System'}; ...
                        vertCsNames(index)];
                end
                index = find(angularUnitsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Angular Units'}; ...
                        angularUnitsNames(index)];
                end
                index = find(linearUnitsCodes == tag);
                if ~isempty(index)
                    nCodesFound = nCodesFound + 1;
                    lasHeaderInfo(36+nCodesFound,1:2) = ...
                        [{'Linear Units'}; ...
                        linearUnitsNames(index)];
                end
            end
        end
    end
end

% Display the header info to the screen in a nice format
disp(['Header information from: ' lasFilename])
for i = 1:length(lasHeaderInfo)
    if isnumeric(lasHeaderInfo{i,2})
        screenOut = [lasHeaderInfo{i,1} ': ' num2str(lasHeaderInfo{i,2})];
    elseif isempty(lasHeaderInfo(i,2))
        screenOut = [lasHeaderInfo{i,1} ': 0'];
    else
        screenOut = [lasHeaderInfo{i,1} ': ' lasHeaderInfo{i,2}];
    end
    disp(screenOut)
end

% Close the LAS file
fclose(fid);



Sable
github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/23136-readlasheader/read_las_header.m
mcbench-benchmarks
github-repos/MATLAB/Sable-mcbench-benchmarks
ba13b2f0296ef49491b95e3f984c7c41fccdb6d8