Index: Damkjer/PointProcessing/SpatialAnalyzer/spanalyze.m =================================================================== --- Damkjer/PointProcessing/SpatialAnalyzer/spanalyze.m (revision 5) +++ Damkjer/PointProcessing/SpatialAnalyzer/spanalyze.m (revision 6) @@ -17,4 +17,6 @@ userParams = parseInputs(varargin{:}); +% disp(userParams); + dimensions=size(QUERIES, 1); elements=size(QUERIES, 2); @@ -33,4 +35,7 @@ step=1000; tic; + +% maxnbrs = 0; + for elem=1:step:elements @@ -38,37 +43,166 @@ % Get the nearest neighbors of elem - if (userParams.neighbors <= 0) - % Perform a fixed radius search - NN=DATABASE.rnn(QUERIES(:,elem:last),... - userParams.radius); + if (userParams.radius > 0 && userParams.neighbors <= 0) + % Perform a fixed radius search + NN=DATABASE.rnn(QUERIES(:,elem:last),... + userParams.radius); + elseif (userParams.radius <= 0 && userParams.neighbors > 0) + % Search unconstrained neighbors + NN=DATABASE.knn(QUERIES(:,elem:last),... + userParams.neighbors); + elseif (userParams.radius > 0 && userParams.neighbors > 0) + % Search constrained to radius + NN=DATABASE.knn(QUERIES(:,elem:last),... + userParams.neighbors,... + 'lim',userParams.radius); + elseif (~isempty(userParams.counts) && userParams.radius <= 0) + % Search unconstrained neighbors + NN=DATABASE.knn(QUERIES(:,elem:last),... + max(userParams.counts)); + elseif (~isempty(userParams.counts) && userParams.radius > 0) + % Search constrained to radius + NN=DATABASE.knn(QUERIES(:,elem:last),... + max(userParams.counts),... + 'lim',userParams.radius); + elseif (~isempty(userParams.steps)) + % Perform a fixed radius search + [NN,DISTS]=DATABASE.rnn(QUERIES(:,elem:last),... + max(userParams.steps)); + end + + if (isempty(userParams.counts) && isempty(userParams.steps)) + [covs,bias,inty]=fastcov(cellfun(@(x) QUERIES(:,x)',NN,'UniformOutput',false)); + elseif (~isempty(userParams.counts)) + tempfeats=zeros(dimensions,size(NN,1)); + tempde=zeros(length(userParams.counts),size(NN,1)); + + tempcovs=cell(length(userParams.counts),1); + tempbias=cell(length(userParams.counts),1); + tempinty=cell(length(userParams.counts),1); + + for c=1:length(userParams.counts) + [tempcovs{c},tempbias{c},tempinty{c}]=fastcov(cellfun(@(x) QUERIES(:,x(1:userParams.counts(c)))',NN,'UniformOutput',false)); + + for nbr=1:size(NN,1) + [~,D]=eig(tempcovs{c}{nbr}); + tempfeats(:,nbr)=sort(diag(D),'descend'); + end + + % Compute dimensional degree + tempdims=(tempfeats(1:dimensions-1,:)-tempfeats(2:dimensions,:))./... + repmat(tempfeats(1,:),dimensions-1,1); + + % Compute isotropy + tempiso=(tempfeats(dimensions,:))./(tempfeats(1,:)); + + % Compute dimensional entropy + tempalpha=[tempdims;tempiso]; + tempde(c,:)=-sum(tempalpha.*log(tempalpha))./log(dimensions); + end + + [~,ind]=min(tempde); + + covs=tempcovs{ind}; + bias=tempbias{ind}; + inty=tempinty{ind}; + elseif (~isempty(userParams.steps)) + tempfeats=zeros(dimensions,size(NN,1)); + tempde=zeros(length(userParams.steps),size(NN,1)); + + tempcovs=cell(length(userParams.steps),1); + tempbias=cell(length(userParams.steps),1); + tempinty=cell(length(userParams.steps),1); + + for c=1:length(userParams.steps) +% disp('userParams.steps(c)'); +% disp(userParams.steps(c)); + +% disp('NN'); +% disp(NN); + +% disp('DISTS'); +% disp(DISTS); + +% disp('DISTS <= userParams.steps(c)'); +% cellfun(@(x) disp(x <= userParams.steps(c))',DISTS,'UniformOutput',false); + + nbrs=cellfun(@(x) sum(x <= userParams.steps(c))',DISTS,'UniformOutput',false); + +% disp('NN(DISTS <= userParams.steps(c))'); +% cellfun(@(x,y) disp(x(y <= userParams.steps(c))'),NN,DISTS,'UniformOutput',false) + +% disp('QUERIES(:,NN(DISTS <= userParams.steps(c)))'); +% cellfun(@(x,y) disp(QUERIES(:,x(y <= userParams.steps(c)))'),NN,DISTS,'UniformOutput',false) + + [tempcovs{c},tempbias{c},tempinty{c}]=fastcov(cellfun(@(x,y) QUERIES(:,x(y <= userParams.steps(c)))',NN,DISTS,'UniformOutput',false)); + +% disp('tempcovs{c}'); +% cellfun(@(x) disp(x), tempcovs{c}); + + for nbr=1:size(NN,1) + if (nbrs{nbr}<5) + continue; + end + + [~,D]=eig(tempcovs{c}{nbr}); + tempfeats(:,nbr)=sort(diag(D),'descend'); + end + +% disp('tempfeats'); +% disp(tempfeats); + + % Compute dimensional degree + % Do not use the following definition based on singular values: + % tempdims=(sqrt(tempfeats(1:dimensions-1,:))-sqrt(tempfeats(2:dimensions,:)))./... + % repmat(sqrt(tempfeats(1,:)),dimensions-1,1); + + tempdims=(tempfeats(1:dimensions-1,:)-tempfeats(2:dimensions,:))./... + repmat(tempfeats(1,:),dimensions-1,1); + + % Compute isotropy + tempiso=(tempfeats(dimensions,:))./(tempfeats(1,:)); + + % Compute dimensional entropy + tempalpha=[tempdims;tempiso]; + tempde(c,:)=-sum(tempalpha.*log(tempalpha))./log(dimensions); + end + +% disp('tempde'); +% disp(tempde); + + [~,ind]=nanmin(tempde); + + covs=tempcovs{ind}; + bias=tempbias{ind}; + inty=tempinty{ind}; + else - if (userParams.radius <= 0) - % Search unconstrained - NN=DATABASE.knn(QUERIES(:,elem:last),... - userParams.neighbors); - else - % Search constrained to radius - NN=DATABASE.knn(QUERIES(:,elem:last),... - userParams.neighbors,... - 'lim',userParams.radius); - end - end - - [covs,bias,inty]=fastcov(cellfun(@(x) QUERIES(:,x)',NN,'UniformOutput',false)); - -% Same solution via SVD appears to be slower via emperical testing + error('What?!'); + end + +% Same solution via SVD appears to be slower based on emperical testing % recentered=fastcenter(cellfun(@(x) QUERIES(:,x)',NN,'UniformOutput',false)); + +% disp('covs'); +% cellfun(@(x) disp(x),covs); for nbr=1:size(NN,1) % skip underconstrained (possible with radius searches) - if (length(NN{nbr})<5) + if (length(NN{nbr})<5 || any(any(isnan(covs{nbr})))) continue; end +% maxnbrs=max(maxnbrs,length(NN{nbr})); + [V,D]=eig(covs{nbr}); [~,index]=min(diag(D)); % norms(:,elem+nbr-1)=abs(V(:,index)) .* sqrt(D(index, index)); norms(:,elem+nbr-1)=abs(V(:,index)); - feats(:,elem+nbr-1)=sort(diag(D),'descend'); + + % Define features as eigenvalues + % feats(:,elem+nbr-1)=sort(diag(D),'descend'); + + % Define features as singular values + feats(:,elem+nbr-1)=sqrt(sort(diag(D),'descend')); biases(elem+nbr-1)=bias{nbr}; @@ -89,4 +223,6 @@ end +% disp(maxnbrs); + if (all(ishghandle(h, 'figure'))) close(h); @@ -120,9 +256,15 @@ % Compute dimensional degree -dims=(feats(1:dimensions-1,:)-feats(2:dimensions,:))./... - repmat(feats(1,:),dimensions-1,1); +dims=zeros(size(feats)); +dims(1:dimensions-1,:)=(feats(1:dimensions-1,:)-feats(2:dimensions,:))./... + repmat(feats(1,:),dimensions-1,1); +dims(dimensions,:)=iso; % Compute Dimensional Embedding [~,embed]=max([dims;iso],[],1); + +% Compute dimensional entropy +alpha=[dims;iso]; +de=-sum(alpha.*log(alpha))./log(dimensions); % Populate feature classes @@ -140,4 +282,5 @@ CLASSES.biases=biases; CLASSES.intensity=ints; +CLASSES.de=de; % st=zeros(size(QUERIES)); @@ -212,5 +355,5 @@ function [userParams] = parseInputs(varargin) -userParams = struct('radius', 0, 'neighbors', 0); +userParams = struct('radius', 0, 'neighbors', 0, 'counts', [], 'steps', []); if length(varargin) == 1 || ~isnumeric(varargin{2}) @@ -237,5 +380,5 @@ end -propertyNames = {'neighbors', 'radius'}; +propertyNames = {'neighbors', 'radius', 'counts', 'steps'}; while ~isempty(varargin) @@ -262,7 +405,7 @@ userParams.neighbors = fix(value); else - error('Damkjer:InvalidCount', ... - ['Number of Neighbors must be a real valued positive '... - 'integer greater or equal to 5']); + error('Damkjer:InvalidCount', ... + ['Number of Neighbors must be a real valued positive '... + 'integer greater or equal to 5']); end case 'radius' @@ -273,7 +416,29 @@ userParams.radius = value; else - error('Damkjer:InvalidPointWeights', ... + error('Damkjer:InvalidRadius', ... 'Radius must be a real valued positive scalar'); end + case 'counts' + if (isvector(value) && ... + isreal(value) && ... + issorted(value) && ... + all(value >= 5)) + userParams.counts = fix(value); + else + error('Damkjer:InvalidCount', ... + ['Counts must be a sorted vector of real valued positive '... + 'integers greater or equal to 5']); + end + case 'steps' + if (isvector(value) && ... + isreal(value) && ... + issorted(value) && ... + all(value > 0)) + userParams.steps = value; + else + error('Damkjer:InvalidSteps', ... + ['Steps must be a sorted vector of real valued positive '... + 'values']); + end end @@ -281,5 +446,23 @@ end -if (userParams.neighbors <= 0 && userParams.radius <= 0) +if (~isempty(userParams.counts) && userParams.neighbors >= 5) + error('Damkjer:MutExOpts', ... + '''neighbors'' and ''counts'' options are mutually exclusive'); +end + +if (~isempty(userParams.steps) && userParams.radius > 0) + error('Damkjer:MutExOpts', ... + '''steps'' and ''radius'' options are mutually exclusive'); +end + +if (~isempty(userParams.counts) && ~isempty(userParams.steps)) + error('Damkjer:MutExOpts', ... + '''steps'' and ''counts'' options are mutually exclusive'); +end + +if (userParams.neighbors <= 0 && ... + userParams.radius <= 0 && ... + isempty(userParams.counts) && ... + isempty(userParams.steps)) userParams.radius = 1; end Index: Damkjer/Util/FileIO/GeoKey/GTModelTypeGeoKey.m =================================================================== --- Damkjer/Util/FileIO/GeoKey/GTModelTypeGeoKey.m (revision 0) +++ Damkjer/Util/FileIO/GeoKey/GTModelTypeGeoKey.m (revision 6) Index: Damkjer/Util/FileIO/GeoKey/GeoKey.m =================================================================== --- Damkjer/Util/FileIO/GeoKey/GeoKey.m (revision 0) +++ Damkjer/Util/FileIO/GeoKey/GeoKey.m (revision 6) @@ -1,2 +1,84 @@ +% GeoKey Abstract class for GeoTIFF Key data structures. +% +% File: +% LasVlr.m +% +% Description: +% This class provides a common ancestor and interface for GeoTIFF Key data +% structures. +% +% Limitations: +% None. +% +% Properties: +% None. +% +% Methods: +% None. +% +% Other m-files required: +% None. +% +% MAT-files required: +% None. +% +% References: +% http://www.remotesensing.org/geotiff/spec/geotiffhome.html +% +% See Also: +% GTModelTypeGeoKey +% GTRasterTypeGeoKey +% GeographicTypeGeoKey +% GeogGeodeticDatumGeoKey +% GeogPrimeMeridianGeoKey +% GeogLinearUnitsGeoKey +% GeogAngularUnitsGeoKey +% GeogEllipsoidGeoKey +% GeogAzimuthUnitsGeoKey +% ProjectedCSTypeGeoKey +% ProjectionGeoKey +% ProjCoordTransGeoKey +% ProjLinearUnitsGeoKey +% VerticalCSTypeGeoKey +% VerticalDatumGeoKey +% VerticalUnitsGeoKey +% GeogLinearUnitSizeGeoKey +% GeogAngularUnitSizeGeoKey +% GeogSemiMajorAxisGeoKey +% GeogSemiMinorAxisGeoKey +% GeogInvFlatteningGeoKey +% GeogPrimeMeridianLongGeoKey +% ProjLinearUnitSizeGeoKey +% ProjStdParallel1GeoKey +% ProjStdParallel2GeoKey +% ProjNatOriginLongGeoKey +% ProjNatOriginLatGeoKey +% ProjFalseEastingGeoKey +% ProjFalseNorthingGeoKey +% ProjFalseOriginLongGeoKey +% ProjFalseOriginLatGeoKey +% ProjFalseOriginEastingGeoKey +% ProjFalseOriginNorthingGeoKey +% ProjCenterLongGeoKey +% ProjCenterLatGeoKey +% ProjCenterEastingGeoKey +% ProjCenterNorthingGeoKey +% ProjScaleAtNatOriginGeoKey +% ProjScaleAtCenterGeoKey +% ProjAzimuthAngleGeoKey +% ProjStraightVertPoleLongGeoKey +% GTCitationGeoKey +% GeogCitationGeoKey +% PCSCitationGeoKey +% VerticalCitationGeoKey +% + +% Software History: +% 2012-AUG-29 K. Damkjer +% Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. +% + classdef GeoKey methods (Abstract = true) Index: Damkjer/Util/FileIO/GeoKey/VerticalUnitsGeoKey.m =================================================================== --- Damkjer/Util/FileIO/GeoKey/VerticalUnitsGeoKey.m (revision 0) +++ Damkjer/Util/FileIO/GeoKey/VerticalUnitsGeoKey.m (revision 6) @@ -1,62 +1,62 @@ classdef (Sealed = true) VerticalUnitsGeoKey < GeoKey - properties - code = 2 - end - - properties (GetAccess = private,Constant = true) - codes = containers.Map(... - {9001,... - 9002,... - 9003,... - 9004,... - 9005,... - 9006,... - 9007,... - 9008,... - 9009,... - 9010,... - 9011,... - 9012,... - 9013,... - 9014,... - 9015},... - {'Meter',... - 'Foot',... - 'Foot (US Survey)',... - 'Foot (Modified American)',... - 'Foot (Clarke)',... - 'Foot (Indian)',... - 'Link',... - 'Link (Benoit)',... - 'Link (Sears)',... - 'Chain (Benoit)',... - 'Chain (Sears)',... - 'Yard (Sears)',... - 'Yard (Indian)',... - 'Fathom',... - 'International Nautical Mile'}); - end - - properties (Dependent) - interpretation - end - - methods - function vlr = VerticalUnitsGeoKey(varargin) - if (nargin == 1) - if (~vlr.codes.isKey(varargin{1})) - error('VerticalUnitsGeoKey:InitError',... - 'Unknown Unit'); - end - vlr.code = varargin{1}; + properties + code = 2 + end + + properties (GetAccess = private,Constant = true) + codes = containers.Map(... + {9001,... + 9002,... + 9003,... + 9004,... + 9005,... + 9006,... + 9007,... + 9008,... + 9009,... + 9010,... + 9011,... + 9012,... + 9013,... + 9014,... + 9015},... + {'Meter',... + 'Foot',... + 'Foot (US Survey)',... + 'Foot (Modified American)',... + 'Foot (Clarke)',... + 'Foot (Indian)',... + 'Link',... + 'Link (Benoit)',... + 'Link (Sears)',... + 'Chain (Benoit)',... + 'Chain (Sears)',... + 'Yard (Sears)',... + 'Yard (Indian)',... + 'Fathom',... + 'International Nautical Mile'}); + end + + properties (Dependent) + interpretation + end + + methods + function vlr = VerticalUnitsGeoKey(varargin) + if (nargin == 1) + if (~vlr.codes.isKey(varargin{1})) + error('VerticalUnitsGeoKey:InitError',... + 'Unknown Unit'); end - end - - function meaning = get.interpretation(vlr) - % Required accessor for dependent property: interpretation - - meaning = vlr.codes(vlr.code); - end - - end + vlr.code = varargin{1}; + end + end + + function meaning = get.interpretation(vlr) + % Required accessor for dependent property: interpretation + + meaning = vlr.codes(vlr.code); + end + + end end Index: Damkjer/Util/FileIO/GeoKey/geoKeyFactory.m =================================================================== --- Damkjer/Util/FileIO/GeoKey/geoKeyFactory.m (revision 0) +++ Damkjer/Util/FileIO/GeoKey/geoKeyFactory.m (revision 6) @@ -1,2 +1,136 @@ +% geoKeyFactory Concrete factory for GeoTIFF key objects +% +% File: +% geoKeyFactory.m +% +% Description: +% This MATLAB function provides an interface similar to a concrete factory +% in the abstract factory design pattern for constructing defined GeoTIFF key +% objects. +% +% Limitations: +% Several of the concrete classes referenced by this factory do not exist. +% This means that an error will be thrown if they key is recognized, but no +% corresponding class exists. This was done to allow me to incrementally +% develop classes for keys as they are encountered. +% +% Synopsis: +% [key] = geoKeyFactory(code, location, size, data) +% +% Inputs: +% vlrHeader - A LAS VLR object containing populated header data +% +% Outputs: +% vlr - A specialized LAS File VLR Specification object +% +% Other m-files required: +% GTModelTypeGeoKey.m +% GTRasterTypeGeoKey.m +% GeographicTypeGeoKey.m +% GeogGeodeticDatumGeoKey.m +% GeogPrimeMeridianGeoKey.m +% GeogLinearUnitsGeoKey.m +% GeogAngularUnitsGeoKey.m +% GeogEllipsoidGeoKey.m +% GeogAzimuthUnitsGeoKey.m +% ProjectedCSTypeGeoKey.m +% ProjectionGeoKey.m +% ProjCoordTransGeoKey.m +% ProjLinearUnitsGeoKey.m +% VerticalCSTypeGeoKey.m +% VerticalDatumGeoKey.m +% VerticalUnitsGeoKey.m +% GeogLinearUnitSizeGeoKey.m +% GeogAngularUnitSizeGeoKey.m +% GeogSemiMajorAxisGeoKey.m +% GeogSemiMinorAxisGeoKey.m +% GeogInvFlatteningGeoKey.m +% GeogPrimeMeridianLongGeoKey.m +% ProjLinearUnitSizeGeoKey.m +% ProjStdParallel1GeoKey.m +% ProjStdParallel2GeoKey.m +% ProjNatOriginLongGeoKey.m +% ProjNatOriginLatGeoKey.m +% ProjFalseEastingGeoKey.m +% ProjFalseNorthingGeoKey.m +% ProjFalseOriginLongGeoKey.m +% ProjFalseOriginLatGeoKey.m +% ProjFalseOriginEastingGeoKey.m +% ProjFalseOriginNorthingGeoKey.m +% ProjCenterLongGeoKey.m +% ProjCenterLatGeoKey.m +% ProjCenterEastingGeoKey.m +% ProjCenterNorthingGeoKey.m +% ProjScaleAtNatOriginGeoKey.m +% ProjScaleAtCenterGeoKey.m +% ProjAzimuthAngleGeoKey.m +% ProjStraightVertPoleLongGeoKey.m +% GTCitationGeoKey.m +% GeogCitationGeoKey.m +% PCSCitationGeoKey.m +% VerticalCitationGeoKey.m +% +% MAT-files required: +% None. +% +% References: +% http://www.remotesensing.org/geotiff/spec/geotiffhome.html +% +% See Also: +% GeoKey +% GTModelTypeGeoKey +% GTRasterTypeGeoKey +% GeographicTypeGeoKey +% GeogGeodeticDatumGeoKey +% GeogPrimeMeridianGeoKey +% GeogLinearUnitsGeoKey +% GeogAngularUnitsGeoKey +% GeogEllipsoidGeoKey +% GeogAzimuthUnitsGeoKey +% ProjectedCSTypeGeoKey +% ProjectionGeoKey +% ProjCoordTransGeoKey +% ProjLinearUnitsGeoKey +% VerticalCSTypeGeoKey +% VerticalDatumGeoKey +% VerticalUnitsGeoKey +% GeogLinearUnitSizeGeoKey +% GeogAngularUnitSizeGeoKey +% GeogSemiMajorAxisGeoKey +% GeogSemiMinorAxisGeoKey +% GeogInvFlatteningGeoKey +% GeogPrimeMeridianLongGeoKey +% ProjLinearUnitSizeGeoKey +% ProjStdParallel1GeoKey +% ProjStdParallel2GeoKey +% ProjNatOriginLongGeoKey +% ProjNatOriginLatGeoKey +% ProjFalseEastingGeoKey +% ProjFalseNorthingGeoKey +% ProjFalseOriginLongGeoKey +% ProjFalseOriginLatGeoKey +% ProjFalseOriginEastingGeoKey +% ProjFalseOriginNorthingGeoKey +% ProjCenterLongGeoKey +% ProjCenterLatGeoKey +% ProjCenterEastingGeoKey +% ProjCenterNorthingGeoKey +% ProjScaleAtNatOriginGeoKey +% ProjScaleAtCenterGeoKey +% ProjAzimuthAngleGeoKey +% ProjStraightVertPoleLongGeoKey +% GTCitationGeoKey +% GeogCitationGeoKey +% PCSCitationGeoKey +% VerticalCitationGeoKey +% + +% Software History: +% 2012-AUG-29 K. Damkjer +% Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. +% + function [key] = geoKeyFactory(code, location, size, data) switch location @@ -36,6 +170,6 @@ key = VerticalUnitsGeoKey(data); otherwise -% error('geoKeyFactory:InvalidGeoNativeCode',... -% 'Unrecognized Code for Native GeoKey'); + error('geoKeyFactory:InvalidGeoNativeCode',... + 'Unrecognized Code for Native GeoKey'); end case 34736 @@ -92,6 +226,6 @@ key = ProjStraightVertPoleLongGeoKey(data); otherwise -% error('geoKeyFactory:InvalidGeoDoubleCode',... -% 'Unrecognized Code for Double GeoKey'); + error('geoKeyFactory:InvalidGeoDoubleCode',... + 'Unrecognized Code for Double GeoKey'); end case 34767 @@ -106,10 +240,10 @@ key = VerticalCitationGeoKey(data,size); otherwise -% error('geoKeyFactory:InvalidGeoAsciiCode',... -% 'Unrecognized Code for Ascii GeoKey'); + error('geoKeyFactory:InvalidGeoAsciiCode',... + 'Unrecognized Code for Ascii GeoKey'); end otherwise -% error('geoKeyFactory:InvalidLocation',... -% 'Unrecognized Data Location for GeoKey'); + error('geoKeyFactory:InvalidLocation',... + 'Unrecognized Data Location for GeoKey'); end end Index: Damkjer/Util/FileIO/LasFile/LasFile.m =================================================================== --- Damkjer/Util/FileIO/LasFile/LasFile.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/LasFile.m (revision 6) @@ -1,613 +1,659 @@ -% LasFile -% -% File: LasFile.m -% -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% file. -% -% Limitations: This class does not currently capture and pass through VLR -% data. There are also currently no mechanisms for setting VLR -% data. This means that it is not possible to geocode LAS -% files with this class. It also means that any geocoding -% information will be lost if data is read in from a geocoded -% file. The output will be invalid unless the geocoding VLRs -% are removed from the LAS file header. +% LasFile ASPRS LAS format data structure +% +% File: +% LasFile.m +% +% Description: +% This MATLAB class represents an ASPRS LAS version 1.2 file. +% +% Limitations: +% There are currently no mechanisms for setting VLR data. This means that it +% is not possible to geocode LAS files with this class. +% +% Properties: +% publicHeader - The LAS File Public Header Block +% vlrs - The LAS File VLRs +% intensity - The Pulse Return Magnitude +% returnNumber - The Pulse Return Number +% numberOfReturns - The Pulse Number of Returns +% scanDirection - The Pulse Scan Direction +% edgeOfFlightLine - Edge of Flight Line Flag +% classification - The Pulse Category Classification +% synthetic - The Synthetic Point Flag +% keyPoint - The Key Point Flag +% withheld - The Withheld Point Flag +% scanAngleRank - The Pulse Scan Angle Rank +% userData - The User Data Field +% pointSourceID - The Point Source ID +% points - The Pulse X, Y, and Z Coordinates +% gpsTime - The Pulse GPS Time +% color - The Point Color +% +% Methods: +% [lasFile] = LasFile(varargin) - Constructor for LasFile objects. +% [lasFile] = loadFrom(location) - Load LAS file from a given location. +% saveTo(location) - Save LAS file to a given location. +% +% Toolbox requirements: +% None. +% +% Script requirements: +% LasPublicHeader.m +% LasVlr.m +% lasVlrFactory.m +% +% Data requirements: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% +% See Also: +% LasPublicHeader +% LasVlr +% lasVlrFactory +% + +% Copyright (C) 2012 Kristian L. Damkjer. % % Software History: % 2012-AUG-29 K. Damkjer % Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % classdef (Sealed = true) LasFile - properties - % The LAS File Public Header Block - % - % The public header block contains information about the LAS file - % that is used to document the details of this file for parsing and - % writing. - publicHeader = LasPublicHeader - - % The LAS File VLRs - % - % The Public Header Block is followed by one or more Variable - % Length Records. Each VLR header is 54 bytes in length with the - % data section dependent on the specific VLR. - vlrs - - % The Pulse Return Magnitude - % - % The intensity value is the integer representation of the pulse - % return magnitude. This value is optional and system specific. - % However, it should always be included if available. - intensity - - % The Pulse Return Number - % - % The Return Number is the pulse return number for a given output - % pulse. A given output laser pulse can have many returns, and they - % must be marked in sequence of return. The first return will have - % a Return Number of one, the second a Return Number of two, and so - % on up to five returns. - returnNumber - - % The Pulse Number of Returns - % - % The Number of Returns is the total number of returns for a given - % pulse. For example, a laser data point may be return two (Return - % Number) within a total of five returns. - numberOfReturns - - % The Pulse Scan Direction - % - % The Scan Direction Flag denotes the direction at which the - % scanner mirror was travelling at the time of the output pulse. A - % bit value of 1 is a positive scan direction, and a bit value of 0 - % is a negative scan direction (where positive scan direction is a - % scan moving from the left side of the in-track direction to the - % right side and negative the opposite). - scanDirection - - % Edge of Flight Line Flag - % - % The Edge of Flight Line data bit has a value of 1 only when the - % point is at the end of a scan. It is the last point on a given - % scan line before it changes direction. - edgeOfFlightLine - - % The Pulse Classification - % - % Classification in LAS 1.0 was essentially user defined and - % optional. Las 1.1 defines a standard set of ASPRS - % classifications. In addition, the field is now mandatory. If a - % point has never been classified, this byte must be set to zero. - % There are no user defined classes since both point format 0 and - % point format 1 supply 8 bits per point for user defined - % operations. - % - % Value | Meaning - % ------+---------------------------------------------------------- - % 0 | Created, never classified - % 1 | Unclassified - % 2 | Ground - % 3 | Low Vegetation - % 4 | Medium Vegetation - % 5 | High Vegetation - % 6 | Building - % 7 | Low Point (noise) - % 8 | Model Key-point (mass point) - % 9 | Water - % 10 | Reserved for ASPRS Definition - % 11 | Reserved for ASPRS Definition - % 12 | Overlap Point - % 13-31 | Reserved for ASPRS Definition - classification - - % The Synthetic Point Flag - % - % If set then this point was created by a technique other than - % LIDAR collection such as digitized from a photogrammetric stereo - % model. - synthetic - - % The Key Point Flag - % - % If set then this point is considered to be a model key-point and - % thus should generally not be withheld in a thinning algorithm. - keyPoint - - % The Withheld Point Flag - % - % If set then this point should not be included in processing - % (synonymous with Deleted). - withheld - - % The Pulse Scan Angle Rank - % - % The Scan Angle Rank is a signed one-byte number with a valid - % range from -90 to +90. The scan angle Rank is the angle (rounded - % to the nearest integer in the absolute value sense) at which the - % laser point was output from the laser system including the roll - % of the aircraft. The scan angle is within 1 degree of accuracy - % from +90 to -90 degrees. The scan angle is an angle based on 0 - % degrees being nadir, and -90 degrees to the left side of the - % aircraft in the direction of flight. - scanAngleRank - - % The User Data Field - % - % This field may be used at the user's discretion. - userData - - % The Point Source ID - % - % This value indicates the file from which this point originated. - % Valid values for this field are 1 to 65,535 inclusive with zero - % being used for a special case discussed below. The numerical - % value corresponds to the File Source ID from which this point - % originated. Zero is reserved as a convenience to system - % implementers. A Point Source ID of zero implies that this point - % originated in this file. This implies that processing software - % should set the Point Source ID equal to the File Source ID of the - % file containing this point at some time during processing. - pointSourceID - end - - properties (Access = private) - % The Pulse X, Y, and Z Coordinates - % - % The X, Y, and Z values are stored as long inegers. The X, Y, and - % Z values are used in conjunction with the scale values and the - % offset values to determine the coordinate for each point as - % described in the Public Header Block section. - pPoints - - % The Pulse GPS Time - % - % The GPS Time is the double floating point time tag value at which - % the point was acquired. It is GPS Week Time if the Global - % Encoding low bit is clear and POSIX Time if the Global Encoding - % low bit is set (see Global Encoding in the Public Header Block - % description). - pGpsTime - - % The Point Color - % - % Red: The Red image channel value associated with this point - % Green: The Green image channel value associated with this point - % Blue: The Blue imagechannel value associated with this point - pColor - end - - properties (Dependent) - % The Pulse X, Y, and Z Coordinates - % - % The X, Y, and Z values are stored as long inegers. The X, Y, and - % Z values are used in conjunction with the scale values and the - % offset values to determine the coordinate for each point as - % described in the Public Header Block section. - points - - % The Pulse GPS Time - % - % The GPS Time is the double floating point time tag value at which - % the point was acquired. It is GPS Week Time if the Global - % Encoding low bit is clear and POSIX Time if the Global Encoding - % low bit is set (see Global Encoding in the Public Header Block - % description). - gpsTime - - % The Point Color - % - % Red: The Red image channel value associated with this point - % Green: The Green image channel value associated with this point - % Blue: The Blue imagechannel value associated with this point - color - end - - methods - function lasFile = LasFile(varargin) - % Constructor for LasFile objects - % - % LasFile objects may be constructed either by using a - % no-argument (default) constructor or by passing a file - % descriptor or file name. + properties + % The LAS File Public Header Block + % + % The public header block contains information about the LAS file + % that is used to document the details of this file for parsing and + % writing. + publicHeader = LasPublicHeader + + % The LAS File VLRs + % + % The Public Header Block is followed by one or more Variable + % Length Records. Each VLR header is 54 bytes in length with the + % data section dependent on the specific VLR. + vlrs + + % The Pulse Return Magnitude + % + % The intensity value is the integer representation of the pulse + % return magnitude. This value is optional and system specific. + % However, it should always be included if available. + intensity + + % The Pulse Return Number + % + % The Return Number is the pulse return number for a given output + % pulse. A given output laser pulse can have many returns, and they + % must be marked in sequence of return. The first return will have + % a Return Number of one, the second a Return Number of two, and so + % on up to five returns. + returnNumber + + % The Pulse Number of Returns + % + % The Number of Returns is the total number of returns for a given + % pulse. For example, a laser data point may be return two (Return + % Number) within a total of five returns. + numberOfReturns + + % The Pulse Scan Direction + % + % The Scan Direction Flag denotes the direction at which the + % scanner mirror was travelling at the time of the output pulse. A + % bit value of 1 is a positive scan direction, and a bit value of 0 + % is a negative scan direction (where positive scan direction is a + % scan moving from the left side of the in-track direction to the + % right side and negative the opposite). + scanDirection + + % Edge of Flight Line Flag + % + % The Edge of Flight Line data bit has a value of 1 only when the + % point is at the end of a scan. It is the last point on a given + % scan line before it changes direction. + edgeOfFlightLine + + % The Pulse Category Classification + % + % Classification in LAS 1.0 was essentially user defined and + % optional. Las 1.1 defines a standard set of ASPRS + % classifications. In addition, the field is now mandatory. If a + % point has never been classified, this byte must be set to zero. + % There are no user defined classes since both point format 0 and + % point format 1 supply 8 bits per point for user defined + % operations. + % + % Value | Meaning + % ------+---------------------------------------------------------- + % 0 | Created, never classified + % 1 | Unclassified + % 2 | Ground + % 3 | Low Vegetation + % 4 | Medium Vegetation + % 5 | High Vegetation + % 6 | Building + % 7 | Low Point (noise) + % 8 | Model Key-point (mass point) + % 9 | Water + % 10 | Reserved for ASPRS Definition + % 11 | Reserved for ASPRS Definition + % 12 | Overlap Point + % 13-31 | Reserved for ASPRS Definition + classification + + % The Synthetic Point Flag + % + % If set then this point was created by a technique other than + % LIDAR collection such as digitized from a photogrammetric stereo + % model. + synthetic + + % The Key Point Flag + % + % If set then this point is considered to be a model key-point and + % thus should generally not be withheld in a thinning algorithm. + keyPoint + + % The Withheld Point Flag + % + % If set then this point should not be included in processing + % (synonymous with Deleted). + withheld + + % The Pulse Scan Angle Rank + % + % The Scan Angle Rank is a signed one-byte number with a valid + % range from -90 to +90. The scan angle Rank is the angle (rounded + % to the nearest integer in the absolute value sense) at which the + % laser point was output from the laser system including the roll + % of the aircraft. The scan angle is within 1 degree of accuracy + % from +90 to -90 degrees. The scan angle is an angle based on 0 + % degrees being nadir, and -90 degrees to the left side of the + % aircraft in the direction of flight. + scanAngleRank + + % The User Data Field + % + % This field may be used at the user's discretion. + userData + + % The Point Source ID + % + % This value indicates the file from which this point originated. + % Valid values for this field are 1 to 65,535 inclusive with zero + % being used for a special case discussed below. The numerical + % value corresponds to the File Source ID from which this point + % originated. Zero is reserved as a convenience to system + % implementers. A Point Source ID of zero implies that this point + % originated in this file. This implies that processing software + % should set the Point Source ID equal to the File Source ID of the + % file containing this point at some time during processing. + pointSourceID + end + + properties (Access = private) + % The Pulse X, Y, and Z Coordinates + % + % The X, Y, and Z values are stored as long inegers. The X, Y, and + % Z values are used in conjunction with the scale values and the + % offset values to determine the coordinate for each point as + % described in the Public Header Block section. + pPoints + + % The Pulse GPS Time + % + % The GPS Time is the double floating point time tag value at which + % the point was acquired. It is GPS Week Time if the Global + % Encoding low bit is clear and POSIX Time if the Global Encoding + % low bit is set (see Global Encoding in the Public Header Block + % description). + pGpsTime + + % The Point Color + % + % Red : The Red image channel value associated with this point + % Green: The Green image channel value associated with this point + % Blue : The Blue imagechannel value associated with this point + pColor + end + + properties (Dependent) + % The Pulse X, Y, and Z Coordinates + % + % The X, Y, and Z values are stored as long integers. The X, Y, and + % Z values are used in conjunction with the scale values and the + % offset values to determine the coordinate for each point as + % described in the Public Header Block section. + points + + % The Pulse GPS Time + % + % The GPS Time is the double floating point time tag value at which + % the point was acquired. It is GPS Week Time if the Global + % Encoding low bit is clear and POSIX Time if the Global Encoding + % low bit is set (see Global Encoding in the Public Header Block + % description). + gpsTime + + % The Point Color + % + % Red : The Red image channel value associated with this point + % Green: The Green image channel value associated with this point + % Blue : The Blue imagechannel value associated with this point + color + end + + methods + function lasFile = LasFile(varargin) + % Constructor for LasFile objects + % + % LasFile objects may be constructed either by using a + % no-argument (default) constructor or by passing a file + % descriptor or file name. + + if (nargin == 0) + lasFile.pPoints=zeros(3,0,'double'); + lasFile.intensity=zeros(1,0,'uint16'); + lasFile.returnNumber=zeros(1,0,'uint8'); + lasFile.numberOfReturns=zeros(1,0,'uint8'); + lasFile.scanDirection=false(1,0); + lasFile.edgeOfFlightLine=false(1,0); + lasFile.classification=zeros(1,0,'uint8'); + lasFile.synthetic=false(1,0); + lasFile.keyPoint=false(1,0); + lasFile.withheld=false(1,0); + lasFile.scanAngleRank=zeros(1,0,'int8'); + lasFile.userData=zeros(1,0,'uint8'); + lasFile.pointSourceID=zeros(1,0,'uint16'); + elseif (nargin == 1) + lasFile = lasFile.loadFrom(varargin{1}); + end + end + + function pointData = get.points(lasFile) + % Required accessor for dependent property: points + pointData = lasFile.pPoints; + end + + function timeData = get.gpsTime(lasFile) + % Required accessor for dependent property: gpsTime + timeData = lasFile.pGpsTime; + end + + function colorData = get.color(lasFile) + % Required accessor for dependent property: color + colorData = lasFile.pColor; + end + + function lasFile = set.points(lasFile, pointData) + % Controlled interface for setting point data on LasFile objects + if (size(pointData,1) ~= 3) + error('LasFile:DimensionMismatch',... + 'Points must have 3 coordinates'); + end + + nPoints = size(pointData,2); + + if (size(lasFile.pPoints,2) > nPoints) + lasFile.intensity=lasFile.intensity(:,1:nPoints); + lasFile.returnNumber=lasFile.returnNumber(:,1:nPoints); + lasFile.numberOfReturns=lasFile.numberOfReturns(:,1:nPoints); + lasFile.scanDirection=lasFile.scanDirection(:,1:nPoints); + lasFile.edgeOfFlightLine=lasFile.edgeOfFlightLine(:,1:nPoints); + lasFile.classification=lasFile.classification(:,1:nPoints); + lasFile.synthetic=lasFile.synthetic(:,1:nPoints); + lasFile.keyPoint=lasFile.keyPoint(:,1:nPoints); + lasFile.withheld=lasFile.withheld(:,1:nPoints); + lasFile.scanAngleRank=lasFile.scanAngleRank(:,1:nPoints); + lasFile.userData=lasFile.userData(:,1:nPoints); + lasFile.pointSourceID=lasFile.pointSourceID(:,1:nPoints); + else + addPoints = nPoints-size(lasFile.pPoints,2); + lasFile.intensity=[lasFile.intensity zeros(1,addPoints,'uint16')]; + lasFile.returnNumber=[lasFile.returnNumber zeros(1,addPoints,'uint8')]; + lasFile.numberOfReturns=[lasFile.numberOfReturns zeros(1,addPoints,'uint8')]; + lasFile.scanDirection=[lasFile.scanDirection false(1,addPoints)]; + lasFile.edgeOfFlightLine=[lasFile.edgeOfFlightLine false(1,addPoints)]; + lasFile.classification=[lasFile.classification zeros(1,addPoints,'uint8')]; + lasFile.synthetic=[lasFile.synthetic false(1,addPoints)]; + lasFile.keyPoint=[lasFile.keyPoint false(1,addPoints)]; + lasFile.withheld=[lasFile.withheld false(1,addPoints)]; + lasFile.scanAngleRank=[lasFile.scanAngleRank zeros(1,addPoints,'int8')]; + lasFile.userData=[lasFile.userData zeros(1,addPoints,'uint8')]; + lasFile.pointSourceID=[lasFile.pointSourceID zeros(1,addPoints,'uint16')]; + end + + lasFile.pPoints = pointData; + + lasFile.publicHeader.nPoints = size(pointData,2); + lasFile.publicHeader.nPointsByReturn(1) = size(pointData,2); + maxX = max(pointData(1,:)); + minX = min(pointData(1,:)); + maxY = max(pointData(2,:)); + minY = min(pointData(2,:)); + maxZ = max(pointData(3,:)); + minZ = min(pointData(3,:)); + lasFile.publicHeader.maxX = maxX; + lasFile.publicHeader.minX = minX; + lasFile.publicHeader.maxY = maxY; + lasFile.publicHeader.minY = minY; + lasFile.publicHeader.maxZ = maxZ; + lasFile.publicHeader.minZ = minZ; + lasFile.publicHeader.xOffset = minX+0.5*(maxX-minX); + lasFile.publicHeader.yOffset = minY+0.5*(maxY-minY); + lasFile.publicHeader.zOffset = minZ+0.5*(maxZ-minZ); + lasFile.publicHeader.xScaleFactor = (maxX-minX)/... + double(intmax('uint32')-1); + lasFile.publicHeader.yScaleFactor = (maxY-minY)/... + double(intmax('uint32')-1); + lasFile.publicHeader.zScaleFactor = (maxZ-minZ)/... + double(intmax('uint32')-1); + end + + function lasFile = set.gpsTime(lasFile, timeData) + % Controlled interface for setting GPS time data on LasFile objects + + if (size(timeData,2) ~= lasFile.publicHeader.nPoints) + error('LasFile:DimensionMismatch',... + 'GPS Time Data must have an entry for each point'); + end + + lasFile.pGpsTime = timeData; + + lasFile.publicHeader.pointDataFormatID = ... + bitor(lasFile.publicHeader.pointDataFormatID,1); + + if (lasFile.publicHeader.pointDataFormatID == 1) + lasFile.publicHeader.pointDataRecordLength = 28; + elseif (lasFile.publicHeader.pointDataFormatID == 3) + lasFile.publicHeader.pointDataRecordLength = 34; + end + end + + function lasFile = set.color(lasFile, colorData) + % Controlled interface for setting color data on LasFile objects. + + if (~strcmp(class(colorData),'uint16')) + error('LasFile:DataTypeMismatch',... + 'Color data must be of type uint16'); + end + + if (size(colorData,1) ~= 3) + error('LasFile:DimensionMismatch',... + 'Color data must have 3 coordinates'); + end + + if (size(colorData,2) ~= lasFile.publicHeader.nPoints) + error('LasFile:DimensionMismatch',... + 'Color data must have an entry for each point'); + end + + lasFile.pColor = colorData; + + lasFile.publicHeader.pointDataFormatID = ... + bitor(lasFile.publicHeader.pointDataFormatID,2); + + if (lasFile.publicHeader.pointDataFormatID == 2) + lasFile.publicHeader.pointDataRecordLength = 26; + elseif (lasFile.publicHeader.pointDataFormatID == 3) + lasFile.publicHeader.pointDataRecordLength = 34; + end + end + + function lasFile = loadFrom(lasFile,location) + % Load the LAS file from the given location. + + if (ischar(location)) + % Assume single argument is file name + [fid, msg]=fopen(location, 'r'); - if (nargin == 0) - lasFile.pPoints=zeros(3,0,'double'); - lasFile.intensity=zeros(1,0,'uint16'); - lasFile.returnNumber=zeros(1,0,'uint8'); - lasFile.numberOfReturns=zeros(1,0,'uint8'); - lasFile.scanDirection=false(1,0); - lasFile.edgeOfFlightLine=false(1,0); - lasFile.classification=zeros(1,0,'uint8'); - lasFile.synthetic=false(1,0); - lasFile.keyPoint=false(1,0); - lasFile.withheld=false(1,0); - lasFile.scanAngleRank=zeros(1,0,'int8'); - lasFile.userData=zeros(1,0,'uint8'); - lasFile.pointSourceID=zeros(1,0,'uint16'); - elseif (nargin == 1) - lasFile = lasFile.loadFrom(varargin{1}); + if (fid < 0) + error('LasFile:FileError',msg); end - end - - function pointData = get.points(lasFile) - % Required accessor for dependent property: points - pointData = lasFile.pPoints; - end - - function timeData = get.gpsTime(lasFile) - % Required accessor for dependent property: gpsTime - timeData = lasFile.pGpsTime; - end - - function colorData = get.color(lasFile) - % Required accessor for dependent property: color - colorData = lasFile.pColor; - end - - function lasFile = set.points(lasFile, pointData) - % Controlled interface for setting point data on LasFile - % objects - if (size(pointData,1) ~= 3) - error('LasFile:DimensionMismatch',... - 'Points must have 3 coordinates'); + elseif (isnumeric(location)) + % Assume single argument is file ID + fid = location; + else + error('LasFile:InitError',... + 'Unknown argument initializer'); + end + + % Parse the LAS Public Header Block + lasFile.publicHeader = LasPublicHeader(fid); + + % Parse the LAS Variable Length Record Headers + if (lasFile.publicHeader.nVLRs > 0) + vlrHeaders = LasVlr(fid); + + my_vlrs = cell(lasFile.publicHeader.nVLRs,1); + + for vlr=1:lasFile.publicHeader.nVLRs + my_vlrs{vlr} = lasVlrFactory(vlrHeaders(vlr)); end - nPoints = size(pointData,2); + lasFile.vlrs = my_vlrs; + end + + % Un-Pack data from uint8 for fast reading. + fseek(fid, lasFile.publicHeader.offsetToPointData, 'bof'); + temp = reshape(fread(fid,... + double(lasFile.publicHeader.pointDataRecordLength)*... + lasFile.publicHeader.nPoints,... + 'uint8=>uint8'),... + lasFile.publicHeader.pointDataRecordLength,... + lasFile.publicHeader.nPoints); + fclose(fid); + + offsets = diag([... + lasFile.publicHeader.xOffset,... + lasFile.publicHeader.yOffset,... + lasFile.publicHeader.zOffset]); + + scales = diag([... + lasFile.publicHeader.xScaleFactor,... + lasFile.publicHeader.yScaleFactor,... + lasFile.publicHeader.zScaleFactor]); + + tempPoints = temp(1:12,:); + tempInt = temp(13:14,:); + tempPtInfo = temp(15,:); + tempPtClass = temp(16,:); + tempSARank = temp(17,:); + tempUData = temp(18,:); + tempSrcID = temp(19:20,:); + + lasFile.pPoints = double(... + reshape(typecast(tempPoints(:),'int32'),... + 3,lasFile.publicHeader.nPoints)); + + lasFile.pPoints = scales*lasFile.pPoints+... + offsets*ones(size(lasFile.pPoints)); + + lasFile.intensity=reshape(typecast(tempInt(:),'uint16'),... + 1,lasFile.publicHeader.nPoints); + + lasFile.returnNumber=bitand(tempPtInfo,7); + lasFile.numberOfReturns=bitand(bitshift(tempPtInfo,-3),7); + lasFile.scanDirection=logical(bitget(tempPtInfo,7)); + lasFile.edgeOfFlightLine=logical(bitget(tempPtInfo,8)); + + lasFile.withheld=logical(bitget(tempPtClass,8)); + lasFile.keyPoint=logical(bitget(tempPtClass,7)); + lasFile.synthetic=logical(bitget(tempPtClass,6)); + lasFile.classification=bitand(tempPtClass,31); + + lasFile.scanAngleRank=reshape(typecast(tempSARank(:),... + 'int8'),1,lasFile.publicHeader.nPoints); + + lasFile.userData=reshape(typecast(tempUData(:),'uint8'),... + 1,lasFile.publicHeader.nPoints); + + lasFile.pointSourceID=reshape(typecast(tempSrcID(:),... + 'uint16'),... + 1,lasFile.publicHeader.nPoints); + + switch lasFile.publicHeader.pointDataFormatID + case 1 + tempGps = temp(21:28,:); + + lasFile.pGpsTime=reshape(typecast(tempGps(:),... + 'double'),... + 1,lasFile.publicHeader.nPoints); + case 2 + tempColor = temp(21:26,:); + + lasFile.pColor = reshape(... + typecast(tempColor(:),'uint16'),... + 3,lasFile.publicHeader.nPoints); + case 3 + tempGps = temp(21:28,:); + tempColor = temp(29:34,:); + + lasFile.pGpsTime=reshape(typecast(tempGps(:),... + 'double'),... + 1,lasFile.publicHeader.nPoints); + + lasFile.pColor = reshape(... + typecast(tempColor(:),'uint16'),... + 3,lasFile.publicHeader.nPoints); + end + + % Update offsets + lasFile.publicHeader.offsetToPointData =... + lasFile.publicHeader.headerSize; + + for record=1:lasFile.publicHeader.nVLRs + lasFile.publicHeader.offsetToPointData = ... + lasFile.publicHeader.offsetToPointData + 54 + ... + lasFile.vlrs{record}.recordLengthAfterHeader; + end + end + + function saveTo(lasFile,location) + % Save the LAS file to the given location. + + if (ischar(location)) + % Assume single argument is file name + [fid, msg]=fopen(location, 'w'); - if (size(lasFile.pPoints,2) > nPoints) - lasFile.intensity=lasFile.intensity(:,1:nPoints); - lasFile.returnNumber=lasFile.returnNumber(:,1:nPoints); - lasFile.numberOfReturns=lasFile.numberOfReturns(:,1:nPoints); - lasFile.scanDirection=lasFile.scanDirection(:,1:nPoints); - lasFile.edgeOfFlightLine=lasFile.edgeOfFlightLine(:,1:nPoints); - lasFile.classification=lasFile.classification(:,1:nPoints); - lasFile.synthetic=lasFile.synthetic(:,1:nPoints); - lasFile.keyPoint=lasFile.keyPoint(:,1:nPoints); - lasFile.withheld=lasFile.withheld(:,1:nPoints); - lasFile.scanAngleRank=lasFile.scanAngleRank(:,1:nPoints); - lasFile.userData=lasFile.userData(:,1:nPoints); - lasFile.pointSourceID=lasFile.pointSourceID(:,1:nPoints); - else - addPoints = nPoints-size(lasFile.pPoints,2); - lasFile.intensity=[lasFile.intensity zeros(1,addPoints,'uint16')]; - lasFile.returnNumber=[lasFile.returnNumber zeros(1,addPoints,'uint8')]; - lasFile.numberOfReturns=[lasFile.numberOfReturns zeros(1,addPoints,'uint8')]; - lasFile.scanDirection=[lasFile.scanDirection false(1,addPoints)]; - lasFile.edgeOfFlightLine=[lasFile.edgeOfFlightLine false(1,addPoints)]; - lasFile.classification=[lasFile.classification zeros(1,addPoints,'uint8')]; - lasFile.synthetic=[lasFile.synthetic false(1,addPoints)]; - lasFile.keyPoint=[lasFile.keyPoint false(1,addPoints)]; - lasFile.withheld=[lasFile.withheld false(1,addPoints)]; - lasFile.scanAngleRank=[lasFile.scanAngleRank zeros(1,addPoints,'int8')]; - lasFile.userData=[lasFile.userData zeros(1,addPoints,'uint8')]; - lasFile.pointSourceID=[lasFile.pointSourceID zeros(1,addPoints,'uint16')]; + if (fid < 0) + error('LasFile:FileError',msg); end + elseif (isnumeric(location)) + % Assume single argument is file ID + fid = location; - lasFile.pPoints = pointData; - - lasFile.publicHeader.nPoints = size(pointData,2); - lasFile.publicHeader.nPointsByReturn(1) = size(pointData,2); - maxX = max(pointData(1,:)); - minX = min(pointData(1,:)); - maxY = max(pointData(2,:)); - minY = min(pointData(2,:)); - maxZ = max(pointData(3,:)); - minZ = min(pointData(3,:)); - lasFile.publicHeader.maxX = maxX; - lasFile.publicHeader.minX = minX; - lasFile.publicHeader.maxY = maxY; - lasFile.publicHeader.minY = minY; - lasFile.publicHeader.maxZ = maxZ; - lasFile.publicHeader.minZ = minZ; - lasFile.publicHeader.xOffset = minX+0.5*(maxX-minX); - lasFile.publicHeader.yOffset = minY+0.5*(maxY-minY); - lasFile.publicHeader.zOffset = minZ+0.5*(maxZ-minZ); - lasFile.publicHeader.xScaleFactor = (maxX-minX)/... - double(intmax('uint32')-1); - lasFile.publicHeader.yScaleFactor = (maxY-minY)/... - double(intmax('uint32')-1); - lasFile.publicHeader.zScaleFactor = (maxZ-minZ)/... - double(intmax('uint32')-1); - end - - function lasFile = set.gpsTime(lasFile, timeData) - % Controlled interface for setting GPS time data on LasFile - % objects - - if (size(timeData,2) ~= lasFile.publicHeader.nPoints) - error('LasFile:DimensionMismatch',... - 'GPS Time Data must have an entry for each point'); - end - - lasFile.pGpsTime = timeData; - - lasFile.publicHeader.pointDataFormatID = ... - bitor(lasFile.publicHeader.pointDataFormatID,1); - - if (lasFile.publicHeader.pointDataFormatID == 1) - lasFile.publicHeader.pointDataRecordLength = 28; - elseif (lasFile.publicHeader.pointDataFormatID == 3) - lasFile.publicHeader.pointDataRecordLength = 34; - end - end - - function lasFile = set.color(lasFile, colorData) - % Controlled interface for setting color data on LasFile - % objects. - - if (~strcmp(class(colorData),'uint16')) - error('LasFile:DataTypeMismatch',... - 'Color data must be of type uint16'); - end - - if (size(colorData,1) ~= 3) - error('LasFile:DimensionMismatch',... - 'Color data must have 3 coordinates'); - end - - if (size(colorData,2) ~= lasFile.publicHeader.nPoints) - error('LasFile:DimensionMismatch',... - 'Color data must have an entry for each point'); - end - - lasFile.pColor = colorData; - - lasFile.publicHeader.pointDataFormatID = ... - bitor(lasFile.publicHeader.pointDataFormatID,2); - - if (lasFile.publicHeader.pointDataFormatID == 2) - lasFile.publicHeader.pointDataRecordLength = 26; - elseif (lasFile.publicHeader.pointDataFormatID == 3) - lasFile.publicHeader.pointDataRecordLength = 34; - end - end - - function lasFile = loadFrom(lasFile,location) - % Load the LAS file from the given location. - - if (ischar(location)) - % Assume single argument is file name - [fid, msg]=fopen(location, 'r'); - - if (fid < 0) - error('LasFile:FileError',msg); - end - elseif (isnumeric(location)) - % Assume single argument is file ID - fid = location; - else - error('LasFile:InitError',... - 'Unknown argument initializer'); - end - - % Parse the LAS Public Header Block - lasFile.publicHeader = LasPublicHeader(fid); - - % Parse the LAS Variable Length Record Headers - if (lasFile.publicHeader.nVLRs > 0) - lasFile.vlrs = LasVlr(fid); - - for vlr=1:lasFile.publicHeader.nVLRs - lasVlrFactory(lasFile.vlrs(vlr)); - end - end - - % Un-Pack data from uint8 for fast reading. - fseek(fid, lasFile.publicHeader.offsetToPointData, 'bof'); - temp = reshape(fread(fid,... - double(lasFile.publicHeader.pointDataRecordLength)*... - lasFile.publicHeader.nPoints,... - 'uint8=>uint8'),... - lasFile.publicHeader.pointDataRecordLength,... - lasFile.publicHeader.nPoints); + else + error('LasFile:InitError',... + 'Unknown argument initializer'); + end + + lasFile.publicHeader.saveTo(fid); + + for record=1:lasFile.publicHeader.nVLRs + lasFile.vlrs{record}.saveTo(fid); + end + + % Pack data to uint8 for fast writing. + offs=diag([... + lasFile.publicHeader.xOffset,... + lasFile.publicHeader.yOffset,... + lasFile.publicHeader.zOffset]); + + invSF=diag([... + 1/lasFile.publicHeader.xScaleFactor,... + 1/lasFile.publicHeader.yScaleFactor,... + 1/lasFile.publicHeader.zScaleFactor]); + + outPoints=int32(invSF*(lasFile.pPoints-... + offs*ones(size(lasFile.pPoints)))); + + uintPoints = reshape(typecast(outPoints(:),'uint8'),... + 12,lasFile.publicHeader.nPoints); + + uintInt = reshape(typecast(lasFile.intensity(:),'uint8'),... + 2,lasFile.publicHeader.nPoints); + + ptInfo=bitand(lasFile.returnNumber,7); + ptInfo=bitor(bitshift(bitand(lasFile.numberOfReturns,7),3),ptInfo); + ptInfo=bitset(ptInfo,7,lasFile.scanDirection); + ptInfo=bitset(ptInfo,8,lasFile.edgeOfFlightLine); + + ptClass=bitand(lasFile.classification,31); + ptClass=bitset(ptClass,6,lasFile.synthetic); + ptClass=bitset(ptClass,7,lasFile.keyPoint); + ptClass=bitset(ptClass,8,lasFile.withheld); + + uintSA = reshape(typecast(lasFile.scanAngleRank(:),'uint8'),... + 1,lasFile.publicHeader.nPoints); + + uintID = reshape(typecast(lasFile.pointSourceID(:),'uint8'),... + 2,lasFile.publicHeader.nPoints); + + switch lasFile.publicHeader.pointDataFormatID + case 0 + temp = [uintPoints;... + uintInt;... + ptInfo;... + ptClass;... + uintSA;... + lasFile.userData;... + uintID]; + case 1 + uintGPS = reshape(typecast(lasFile.pGpsTime(:),'uint8'),... + 8,lasFile.publicHeader.nPoints); + + temp = [uintPoints;... + uintInt;... + ptInfo;... + ptClass;... + uintSA;... + lasFile.userData;... + uintID;... + uintGPS]; + case 2 + uintColor = reshape(typecast(lasFile.pColor(:),'uint8'),... + 6,lasFile.publicHeader.nPoints); + + temp = [uintPoints;... + uintInt;... + ptInfo;... + ptClass;... + uintSA;... + lasFile.userData;... + uintID;... + uintColor]; + case 3 + uintGPS = reshape(typecast(lasFile.pGpsTime(:),'uint8'),... + 8,lasFile.publicHeader.nPoints); + + uintColor = reshape(typecast(lasFile.pColor(:),'uint8'),... + 6,lasFile.publicHeader.nPoints); + + temp = [uintPoints;... + uintInt;... + ptInfo;... + ptClass;... + uintSA;... + lasFile.userData;... + uintID;... + uintGPS;... + uintColor]; + end + + fseek(fid, lasFile.publicHeader.offsetToPointData, 'bof'); + fwrite(fid,temp); + + if (ischar(location)) + % Close the file if it was opened here. fclose(fid); - - offsets = diag([... - lasFile.publicHeader.xOffset,... - lasFile.publicHeader.yOffset,... - lasFile.publicHeader.zOffset]); - - scales = diag([... - lasFile.publicHeader.xScaleFactor,... - lasFile.publicHeader.yScaleFactor,... - lasFile.publicHeader.zScaleFactor]); - - tempPoints = temp(1:12,:); - tempInt = temp(13:14,:); - tempPtInfo = temp(15,:); - tempPtClass = temp(16,:); - tempSARank = temp(17,:); - tempUData = temp(18,:); - tempSrcID = temp(19:20,:); - - lasFile.pPoints = double(... - reshape(typecast(tempPoints(:),'int32'),... - 3,lasFile.publicHeader.nPoints)); - - lasFile.pPoints = scales*lasFile.pPoints+... - offsets*ones(size(lasFile.pPoints)); - - lasFile.intensity=reshape(typecast(tempInt(:),'uint16'),... - 1,lasFile.publicHeader.nPoints); - - lasFile.returnNumber=bitand(tempPtInfo,7); - lasFile.numberOfReturns=bitand(bitshift(tempPtInfo,-3),7); - lasFile.scanDirection=logical(bitget(tempPtInfo,7)); - lasFile.edgeOfFlightLine=logical(bitget(tempPtInfo,8)); - - lasFile.withheld=logical(bitget(tempPtClass,8)); - lasFile.keyPoint=logical(bitget(tempPtClass,7)); - lasFile.synthetic=logical(bitget(tempPtClass,6)); - lasFile.classification=bitand(tempPtClass,31); - - lasFile.scanAngleRank=reshape(typecast(tempSARank(:),... - 'int8'),1,lasFile.publicHeader.nPoints); - - lasFile.userData=reshape(typecast(tempUData(:),'uint8'),... - 1,lasFile.publicHeader.nPoints); - - lasFile.pointSourceID=reshape(typecast(tempSrcID(:),... - 'uint16'),... - 1,lasFile.publicHeader.nPoints); - - switch lasFile.publicHeader.pointDataFormatID - case 1 - tempGps = temp(21:28,:); - - lasFile.pGpsTime=reshape(typecast(tempGps(:),... - 'double'),... - 1,lasFile.publicHeader.nPoints); - case 2 - tempColor = temp(21:26,:); - - lasFile.pColor = reshape(... - typecast(tempColor(:),'uint16'),... - 3,lasFile.publicHeader.nPoints); - case 3 - tempGps = temp(21:28,:); - tempColor = temp(29:34,:); - - lasFile.pGpsTime=reshape(typecast(tempGps(:),... - 'double'),... - 1,lasFile.publicHeader.nPoints); - - lasFile.pColor = reshape(... - typecast(tempColor(:),'uint16'),... - 3,lasFile.publicHeader.nPoints); - end - - % Update offsets - lasFile.publicHeader.offsetToPointData =... - lasFile.publicHeader.headerSize; - - for record=1:lasFile.publicHeader.nVLRs - lasFile.publicHeader.offsetToPointData = ... - lasFile.publicHeader.offsetToPointData + 54 + ... - lasFile.vlrs(record).recordLengthAfterHeader; - end - end - - function saveTo(lasFile,location) - % Save the LAS file to the given location. - - if (ischar(location)) - % Assume single argument is file name - [fid, msg]=fopen(location, 'w'); - - if (fid < 0) - error('LasFile:FileError',msg); - end - elseif (isnumeric(location)) - % Assume single argument is file ID - fid = location; - - else - error('LasFile:InitError',... - 'Unknown argument initializer'); - end - - lasFile.publicHeader.saveTo(fid); - - for record=1:lasFile.publicHeader.nVLRs - lasFile.vlrs(record).saveTo(fid); - end - - % Pack data to uint8 for fast writing. - offs=diag([... - lasFile.publicHeader.xOffset,... - lasFile.publicHeader.yOffset,... - lasFile.publicHeader.zOffset]); - - invSF=diag([... - 1/lasFile.publicHeader.xScaleFactor,... - 1/lasFile.publicHeader.yScaleFactor,... - 1/lasFile.publicHeader.zScaleFactor]); - - outPoints=int32(invSF*(lasFile.pPoints-... - offs*ones(size(lasFile.pPoints)))); - - uintPoints = reshape(typecast(outPoints(:),'uint8'),... - 12,lasFile.publicHeader.nPoints); - - uintInt = reshape(typecast(lasFile.intensity(:),'uint8'),... - 2,lasFile.publicHeader.nPoints); - - ptInfo=bitand(lasFile.returnNumber,7); - ptInfo=bitor(bitshift(bitand(lasFile.numberOfReturns,7),3),ptInfo); - ptInfo=bitset(ptInfo,7,lasFile.scanDirection); - ptInfo=bitset(ptInfo,8,lasFile.edgeOfFlightLine); - - ptClass=bitand(lasFile.classification,31); - ptClass=bitset(ptClass,6,lasFile.synthetic); - ptClass=bitset(ptClass,7,lasFile.keyPoint); - ptClass=bitset(ptClass,8,lasFile.withheld); - - uintSA = reshape(typecast(lasFile.scanAngleRank(:),'uint8'),... - 1,lasFile.publicHeader.nPoints); - - uintID = reshape(typecast(lasFile.pointSourceID(:),'uint8'),... - 2,lasFile.publicHeader.nPoints); - - switch lasFile.publicHeader.pointDataFormatID - case 0 - temp = [uintPoints;... - uintInt;... - ptInfo;... - ptClass;... - uintSA;... - lasFile.userData;... - uintID]; - case 1 - uintGPS = reshape(typecast(lasFile.pGpsTime(:),'uint8'),... - 8,lasFile.publicHeader.nPoints); - - temp = [uintPoints;... - uintInt;... - ptInfo;... - ptClass;... - uintSA;... - lasFile.userData;... - uintID;... - uintGPS]; - case 2 - uintColor = reshape(typecast(lasFile.pColor(:),'uint8'),... - 6,lasFile.publicHeader.nPoints); - - temp = [uintPoints;... - uintInt;... - ptInfo;... - ptClass;... - uintSA;... - lasFile.userData;... - uintID;... - uintColor]; - case 3 - uintGPS = reshape(typecast(lasFile.pGpsTime(:),'uint8'),... - 8,lasFile.publicHeader.nPoints); - - uintColor = reshape(typecast(lasFile.pColor(:),'uint8'),... - 6,lasFile.publicHeader.nPoints); - - temp = [uintPoints;... - uintInt;... - ptInfo;... - ptClass;... - uintSA;... - lasFile.userData;... - uintID;... - uintGPS;... - uintColor]; - end - - fseek(fid, lasFile.publicHeader.offsetToPointData, 'bof'); - fwrite(fid,temp); - - if (ischar(location)) - % Close the file if it was opened here. - fclose(fid); - end - end - end + end + end + end end Index: Damkjer/Util/FileIO/LasFile/LasPublicHeader.m =================================================================== --- Damkjer/Util/FileIO/LasFile/LasPublicHeader.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/LasPublicHeader.m (revision 6) @@ -1,580 +1,643 @@ -% LasPublicHeader -% -% File: LasPublicHeader.m -% -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% file public header block. -% -% Software History: -% 2012-AUG-29 K. Damkjer -% Initial Coding. +% LasPublicHeader ASPRS LAS format public header block +% +% File: +% LasPublicHeader.m +% +% Description: +% This MATLAB class represents an ASPRS LAS version 1.2 file public header +% block. +% +% Limitations: +% None. +% +% Properties: +% fileSignature - The File Signature +% fileSourceID - The File Source ID +% globalEncoding - The File Global Encoding +% projectID_GUID1 - The Globally Unique (Project) Identifier - Part 1 +% projectID_GUID2 - The Globally Unique (Project) Identifier - Part 2 +% projectID_GUID3 - The Globally Unique (Project) Identifier - Part 3 +% projectID_GUID4 - The Globally Unique (Project) Identifier - Part 4 +% systemIdentifier - The System Identifier +% generatingSoftware - The Generating Software +% headerSize - The Header Size in Bytes +% offsetToPointData - The Offset to Point Data +% nVLRs - The Number of Variable Length Records +% pointDataFormatID - The Point Data Format ID +% pointDataRecordLength - The Point Data Record Length +% nPoints - The Number of Point Records +% nPointsByReturn - The Number of Point Records by Return +% xScaleFactor - The X Coordinate Scale Factor +% yScaleFactor - The Y Coordinate Scale Factor +% zScaleFactor - The Z Coordinate Scale Factor +% xOffset - The X Coordinate Offset +% yOffset - The Y Coordinate Offset +% zOffset - The Z Coordinate Offset +% maxX - The Maximum X Coordinate +% minX - The Minimum X Coordinate +% maxY - The Maximum Y Coordinate +% minY - The Minimum Y Coordinate +% maxZ - The Maximum Z Coordinate +% minZ - The Minimum Z Coordinate +% extraData - The Extra Data +% fileCreationDate - The File Creation Date +% fileCreationDayOfYear - The File Creation Day of Year +% fileCreationYear - The File Creation Year +% version - The Version Number +% +% Methods: +% [header] = LasPublicHeader(varargin) - Constructor for LasPublicHeader +% objects. +% [header] = loadFrom(location) - Load LAS file public header block +% from a given location. +% saveTo(location) - Save LAS file public header block +% to a given location. +% +% Other m-files required: +% None. +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% +% See Also: +% LasFile +% + +% Software History: +% 2012-AUG-29 K. Damkjer +% Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % classdef (Sealed = true) LasPublicHeader - properties (SetAccess = private) - % The File Signature - % The file signature must contain the four characters "LASF", and - % it is required by the LAS specification. These four characters - % can be checked by user software as a quick look initial - % determination of file type. - fileSignature = 'LASF' - end - - properties - % The File Source ID - % - % This field should be set to a value between 1 and 65535, - % inclusive. A value of zero (0) is interpreted to mean that an ID - % has not been assigned. In this case, processing software is free - % to assign any valid number. Note that this scheme allows a LIDAR - % project to contain up to 65535 unique sources. A source can be - % considered an original flight line or it can be the result of - % merge and/or extract operations. - fileSourceID = uint16(0) - - % The File Global Encoding - % - % This is a bit field used to indicate certain global properties - % about the file. In LAS 1.2 (the version in which this field is - % introduced), only the low bit is defined (this is the bit, that - % if set, would have the unsigned integer yield a value of 1). This - % bit field is defined as: - % - % Bits | Field Name | Description - % -----+---------------+----------------------------------------- - % 0 | GPS Time Type | The meaning of GPS Time in the Point - % | | Records - % | | 0(not set) -> GPS Time in the point - % | | record fields is GPS Week Time (the same - % | | as previous versions of LAS) - % | | 1 (set) -> GPS Time is standard GPS Time - % | | (satellite GPS Time) minus 1 x 10^9. The - % | | offset moves the time back to near zero - % | | to improve floating point resolution. - % -----+---------------+----------------------------------------- - % 1:15 | Reserved | Must be set to zero. - globalEncoding = uint16(0) - - % The Globally Unique (Project) Identifier - Part 1 - % - % The four fields that comprise a complete Globally Unique - % Identifier (GUID) are now reserved for use as a Project - % Identifier (Project ID). The field remains optional. The time of - % assignment of the Project ID is at the discretion of processing - % software. The Project ID should be the same for all files that - % are associated with a unique project. By assigning a Project ID - % and using a File Source ID (defined above) every file within a - % project and every point within a file can be uniquely identified, - % globally. - projectID_GUID1 = uint32(0) - - % The Globally Unique (Project) Identifier - Part 2 - % - % The four fields that comprise a complete Globally Unique - % Identifier (GUID) are now reserved for use as a Project - % Identifier (Project ID). The field remains optional. The time of - % assignment of the Project ID is at the discretion of processing - % software. The Project ID should be the same for all files that - % are associated with a unique project. By assigning a Project ID - % and using a File Source ID (defined above) every file within a - % project and every point within a file can be uniquely identified, - % globally. - projectID_GUID2 = uint16(0) - - % The Globally Unique (Project) Identifier - Part 3 - % - % The four fields that comprise a complete Globally Unique - % Identifier (GUID) are now reserved for use as a Project - % Identifier (Project ID). The field remains optional. The time of - % assignment of the Project ID is at the discretion of processing - % software. The Project ID should be the same for all files that - % are associated with a unique project. By assigning a Project ID - % and using a File Source ID (defined above) every file within a - % project and every point within a file can be uniquely identified, - % globally. - projectID_GUID3 = uint16(0) - - % The Globally Unique (Project) Identifier - Part 4 - % - % The four fields that comprise a complete Globally Unique - % Identifier (GUID) are now reserved for use as a Project - % Identifier (Project ID). The field remains optional. The time of - % assignment of the Project ID is at the discretion of processing - % software. The Project ID should be the same for all files that - % are associated with a unique project. By assigning a Project ID - % and using a File Source ID (defined above) every file within a - % project and every point within a file can be uniquely identified, - % globally. - projectID_GUID4 = '' - - % The System Identifier - % - % The version 1.0 specification assumes that LAS files are - % exclusively generated as a result of collection by a hardware - % sensor. Version 1.1 recognizes that files often result from - % extraction, merging or modifying existing data files. Thus System - % ID becomes: - % - % Generating Agent | System ID - % ---------------------------------------+----------------------- - % Hardware System | String identifying - % | hardware (e.g. "ALTM - % | 1210" or "ALS50" - % ---------------------------------------+----------------------- - % Merge of one or more files | "MERGE" - % ---------------------------------------+----------------------- - % Modification of a single file | "MODIFICATION" - % ---------------------------------------+----------------------- - % Extraction from one or more files | "EXTRACTION" - % ---------------------------------------+----------------------- - % Reprojection, rescaling, warping, etc. | "TRANSFORMATION" - % ---------------------------------------+----------------------- - % Some other operation | "OTHER" or a string up - % | to 32 characters - % | identifying the - % | operation - systemIdentifier = 'OTHER' - - % The Generating Software - % - % This information is ASCII data describing the generating software - % itself. This field provides a mechanism for specifying which - % generating software package and version was used during LAS file - % creation (e.g. "TerraScan V-10.8", "REALM V-4.2" and etc.). If - % the character data is less than 32 characters, the remaining data - % must be null. - generatingSoftware = ['MATLAB ' version] - - % The Header Size - % - % The size, in bytes, of the Public Header Block itself. In the - % event that the header is extended by a software application - % through the addition of data at the end of the header, the Header - % Size field must be updated with the new header size. Extension of - % the Public Header Block is discouraged; the Variable Length - % Records should be used whenever possible to add custom header - % data. In the event a generating software package adds data to the - % Public Header Block, this data must be placed at the end of the - % structure and the Header Size must be updated to reflect the new - % size. - headerSize = uint16(227) - - % The Offset to Point Data - % - % The actual number of bytes from the beginning of the file to the - % first field of the first point record data field. This data - % offset must be updated if any software adds data from the Public - % Header Block or adds/removes data to/from the Variable Length - % Records. - offsetToPointData = uint32(227) - - % The Number of Variable Length Records - % This field contains the current number of Variable Length - % Records. This number must be updated if the number of Variable - % Length Records changes at any time. - nVLRs = uint32(0) - - % The Point Data Format ID - % - % The point data format ID corresponds to the point data record - % format type. LAS 1.2 defines types 0, 1, 2 and 3. - pointDataFormatID = uint8(0) - - % The Point Data Record Length - % - % The size, in bytes, of the Point Data Record. - pointDataRecordLength = uint16(20) - - % The Number of Point Records - % - % This field contains the total number of point records within the - % file. - nPoints = uint32(0) - - % The Number of Points by Return - % This field contains an array of the total point records per - % return. The first unsigned long value will be the total number of - % records from the first return, and the second contains the total - % number for return two, and so forth up to five returns. - nPointsByReturn = zeros(5, 1, 'uint32') - - % The X Coordinate Scale Factor - % - % The scale factor fields contain a double floating point value - % that is uysed to scale the corresponding X, Y, and Z long values - % within the point records. The corresponding X, Y, and Z scale - % factor must be multiplied by the X, Y, or Z point record value to - % get the actual X, Y, or Z coordinate. For example, if the X, Y, - % and Z coordinates are intended to have two decimal point values, - % then each scale factor will contain the number 0.01. - xScaleFactor = double(0) + properties (SetAccess = private) + % The File Signature + % + % The file signature must contain the four characters "LASF", and + % it is required by the LAS specification. These four characters + % can be checked by user software as a quick look initial + % determination of file type. + fileSignature = 'LASF' + end + + properties + % The File Source ID + % + % This field should be set to a value between 1 and 65535, + % inclusive. A value of zero (0) is interpreted to mean that an ID + % has not been assigned. In this case, processing software is free + % to assign any valid number. Note that this scheme allows a LIDAR + % project to contain up to 65535 unique sources. A source can be + % considered an original flight line or it can be the result of + % merge and/or extract operations. + fileSourceID = uint16(0) + + % The File Global Encoding + % + % This is a bit field used to indicate certain global properties + % about the file. In LAS 1.2 (the version in which this field is + % introduced), only the low bit is defined (this is the bit, that + % if set, would have the unsigned integer yield a value of 1). This + % bit field is defined as: + % + % Bits | Field Name | Description + % -----+---------------+----------------------------------------- + % 0 | GPS Time Type | The meaning of GPS Time in the Point + % | | Records + % | | 0(not set) -> GPS Time in the point + % | | record fields is GPS Week Time (the same + % | | as previous versions of LAS) + % | | 1 (set) -> GPS Time is standard GPS Time + % | | (satellite GPS Time) minus 1 x 10^9. The + % | | offset moves the time back to near zero + % | | to improve floating point resolution. + % -----+---------------+----------------------------------------- + % 1:15 | Reserved | Must be set to zero. + globalEncoding = uint16(0) + + % The Globally Unique (Project) Identifier - Part 1 + % + % The four fields that comprise a complete Globally Unique + % Identifier (GUID) are now reserved for use as a Project + % Identifier (Project ID). The field remains optional. The time of + % assignment of the Project ID is at the discretion of processing + % software. The Project ID should be the same for all files that + % are associated with a unique project. By assigning a Project ID + % and using a File Source ID (defined above) every file within a + % project and every point within a file can be uniquely identified, + % globally. + projectID_GUID1 = uint32(0) + + % The Globally Unique (Project) Identifier - Part 2 + % + % The four fields that comprise a complete Globally Unique + % Identifier (GUID) are now reserved for use as a Project + % Identifier (Project ID). The field remains optional. The time of + % assignment of the Project ID is at the discretion of processing + % software. The Project ID should be the same for all files that + % are associated with a unique project. By assigning a Project ID + % and using a File Source ID (defined above) every file within a + % project and every point within a file can be uniquely identified, + % globally. + projectID_GUID2 = uint16(0) + + % The Globally Unique (Project) Identifier - Part 3 + % + % The four fields that comprise a complete Globally Unique + % Identifier (GUID) are now reserved for use as a Project + % Identifier (Project ID). The field remains optional. The time of + % assignment of the Project ID is at the discretion of processing + % software. The Project ID should be the same for all files that + % are associated with a unique project. By assigning a Project ID + % and using a File Source ID (defined above) every file within a + % project and every point within a file can be uniquely identified, + % globally. + projectID_GUID3 = uint16(0) + + % The Globally Unique (Project) Identifier - Part 4 + % + % The four fields that comprise a complete Globally Unique + % Identifier (GUID) are now reserved for use as a Project + % Identifier (Project ID). The field remains optional. The time of + % assignment of the Project ID is at the discretion of processing + % software. The Project ID should be the same for all files that + % are associated with a unique project. By assigning a Project ID + % and using a File Source ID (defined above) every file within a + % project and every point within a file can be uniquely identified, + % globally. + projectID_GUID4 = '' + + % The System Identifier + % + % The version 1.0 specification assumes that LAS files are + % exclusively generated as a result of collection by a hardware + % sensor. Version 1.1 recognizes that files often result from + % extraction, merging or modifying existing data files. Thus System + % ID becomes: + % + % Generating Agent | System ID + % ---------------------------------------+----------------------- + % Hardware System | String identifying + % | hardware (e.g. "ALTM + % | 1210" or "ALS50" + % ---------------------------------------+----------------------- + % Merge of one or more files | "MERGE" + % ---------------------------------------+----------------------- + % Modification of a single file | "MODIFICATION" + % ---------------------------------------+----------------------- + % Extraction from one or more files | "EXTRACTION" + % ---------------------------------------+----------------------- + % Reprojection, rescaling, warping, etc. | "TRANSFORMATION" + % ---------------------------------------+----------------------- + % Some other operation | "OTHER" or a string up + % | to 32 characters + % | identifying the + % | operation + systemIdentifier = 'OTHER' + + % The Generating Software + % + % This information is ASCII data describing the generating software + % itself. This field provides a mechanism for specifying which + % generating software package and version was used during LAS file + % creation (e.g. "TerraScan V-10.8", "REALM V-4.2" and etc.). If + % the character data is less than 32 characters, the remaining data + % must be null. + generatingSoftware = ['MATLAB ' version] + + % The Header Size + % + % The size, in bytes, of the Public Header Block itself. In the + % event that the header is extended by a software application + % through the addition of data at the end of the header, the Header + % Size field must be updated with the new header size. Extension of + % the Public Header Block is discouraged; the Variable Length + % Records should be used whenever possible to add custom header + % data. In the event a generating software package adds data to the + % Public Header Block, this data must be placed at the end of the + % structure and the Header Size must be updated to reflect the new + % size. + headerSize = uint16(227) + + % The Offset to Point Data + % + % The actual number of bytes from the beginning of the file to the + % first field of the first point record data field. This data + % offset must be updated if any software adds data from the Public + % Header Block or adds/removes data to/from the Variable Length + % Records. + offsetToPointData = uint32(227) + + % The Number of Variable Length Records + % This field contains the current number of Variable Length + % Records. This number must be updated if the number of Variable + % Length Records changes at any time. + nVLRs = uint32(0) + + % The Point Data Format ID + % + % The point data format ID corresponds to the point data record + % format type. LAS 1.2 defines types 0, 1, 2 and 3. + pointDataFormatID = uint8(0) + + % The Point Data Record Length + % + % The size, in bytes, of the Point Data Record. + pointDataRecordLength = uint16(20) + + % The Number of Point Records + % + % This field contains the total number of point records within the + % file. + nPoints = uint32(0) + + % The Number of Points by Return + % + % This field contains an array of the total point records per + % return. The first unsigned long value will be the total number of + % records from the first return, and the second contains the total + % number for return two, and so forth up to five returns. + nPointsByReturn = zeros(5, 1, 'uint32') + + % The X Coordinate Scale Factor + % + % The scale factor fields contain a double floating point value + % that is uysed to scale the corresponding X, Y, and Z long values + % within the point records. The corresponding X, Y, and Z scale + % factor must be multiplied by the X, Y, or Z point record value to + % get the actual X, Y, or Z coordinate. For example, if the X, Y, + % and Z coordinates are intended to have two decimal point values, + % then each scale factor will contain the number 0.01. + xScaleFactor = double(0) + + % The Y Coordinate Scale Factor + % + % The scale factor fields contain a double floating point value + % that is uysed to scale the corresponding X, Y, and Z long values + % within the point records. The corresponding X, Y, and Z scale + % factor must be multiplied by the X, Y, or Z point record value to + % get the actual X, Y, or Z coordinate. For example, if the X, Y, + % and Z coordinates are intended to have two decimal point values, + % then each scale factor will contain the number 0.01. + yScaleFactor = double(0) + + % The Z Coordinate Scale Factor + % + % The scale factor fields contain a double floating point value + % that is uysed to scale the corresponding X, Y, and Z long values + % within the point records. The corresponding X, Y, and Z scale + % factor must be multiplied by the X, Y, or Z point record value to + % get the actual X, Y, or Z coordinate. For example, if the X, Y, + % and Z coordinates are intended to have two decimal point values, + % then each scale factor will contain the number 0.01. + zScaleFactor = double(0) + + % The X Coordinate Offset + % + % The offset fields should be used to set the overall offset for + % the point records. In general these numbers will be zero, but for + % certain cases the resolution of the point data may not be large + % enough for a given projection system. However, it should always + % be assumed that these numbers are used. So to scale a given X + % from the point record, take the point record X multiplied by the + % X scale factor, and then add the X offset. + % + % X_coordinate = (X_record * X_scale) + X_offset + xOffset = double(0) + + % The Y Coordinate Offset + % + % The offset fields should be used to set the overall offset for + % the point records. In general these numbers will be zero, but for + % certain cases the resolution of the point data may not be large + % enough for a given projection system. However, it should always + % be assumed that these numbers are used. So to scale a given Y + % from the point record, take the point record Y multiplied by the + % Y scale factor, and then add the Y offset. + % + % Y_coordinate = (Y_record * Y_scale) + Y_offset + yOffset = double(0) + + % The Z Coordinate Offset + % + % The offset fields should be used to set the overall offset for + % the point records. In general these numbers will be zero, but for + % certain cases the resolution of the point data may not be large + % enough for a given projection system. However, it should always + % be assumed that these numbers are used. So to scale a given Z + % from the point record, take the point record Z multiplied by the + % Z scale factor, and then add the Z offset. + % + % Z_coordinate = (Z_record * Z_scale) + Z_offset + zOffset = double(0) + + % The Maximum X Coordinate + % + % The max and min data fields are the actual unscaled extents of + % the LAS point file data, specified in the coordinate system of + % the LAS data. + maxX = double(0) + + % The Minimum X Coordinate + % + % The max and min data fields are the actual unscaled extents of + % the LAS point file data, specified in the coordinate system of + % the LAS data. + minX = double(0) + + % The Maximum Y Coordinate + % + % The max and min data fields are the actual unscaled extents of + % the LAS point file data, specified in the coordinate system of + % the LAS data. + maxY = double(0) + + % The Minimum Y Coordinate + % + % The max and min data fields are the actual unscaled extents of + % the LAS point file data, specified in the coordinate system of + % the LAS data. + minY = double(0) + + % The Maximum Z Coordinate + % + % The max and min data fields are the actual unscaled extents of + % the LAS point file data, specified in the coordinate system of + % the LAS data. + maxZ = double(0) + + % The Minimum Z Coordinate + % + % The max and min data fields are the actual unscaled extents of + % the LAS point file data, specified in the coordinate system of + % the LAS data. + minZ = double(0) + + % The Extra Data + % + % The public header may be extended, though it is highly + % discouraged. If it is extended, we have no way of knowing how to + % interpret the extra data. We can simply pass it through, though. + extraData = zeros(0,'uint8'); + end + + properties (Access = private) + % The File Creation Date + % + % The date on which the file was created expressed as a serial date + % number. + pFileCreationDate + + % The Major Version Number + % + % The version number consists of a major and minor field. The major + % and minor fields combine to form the number that indicates the + % format number of the current specification itself. For example, + % specification number 1.2 (this version) would contain 1 in the + % major field and 2 in the minor field. + pVersionMajor = uint8(1) + + % The Minor Version Number + % + % The version number consists of a major and minor field. The major + % and minor fields combine to form the number that indicates the + % format number of the current specification itself. For example, + % specification number 1.2 (this version) would contain 1 in the + % major field and 2 in the minor field. + pVersionMinor = uint8(2) + end + + properties (Dependent) + % The File Creation Date + % + % The date on which the file was created expressed as a date + % string. + fileCreationDate + + % The File Creation Day of Year + % + % The day, expressed as an unsigned short, on which this file was + % created. Day is computed as the Greenwhich Mean Time (GMT) day. + % January 1 is considered day 1. + fileCreationDayOfYear + + % The File Creation Year + % + % The year, expressed as a four digit number, in which the file was + % created. + fileCreationYear + + % The Version Number + % + % The version number consists of a major and minor field. The major + % and minor fields combine to form the number that indicates the + % format number of the current specification itself. For example, + % specification number 1.2 (this version) would contain 1 in the + % major field and 2 in the minor field. + version + end + + methods + function header = LasPublicHeader(varargin) + % Constructor for LasFile public header block objects + % + % LasPublicHeader objects may be constructed either by using a + % no-argument (default) constructor or by passing a file + % descriptor or file name. + switch (nargin) + case 0 + header.pFileCreationDate = floor(datenum(1970,1,1)+... + java.lang.System.currentTimeMillis/... + (1e3*86400)); + case 1 + header = header.loadFrom(varargin{1}); + otherwise + error('LasPublicHeader:UnexpectedInputs',... + 'Unexpected number of inputs encountered.'); + end + end + + function theDate = get.fileCreationDate(header) + theDate = datestr(header.pFileCreationDate); + end + + function theDOY = get.fileCreationDayOfYear(header) + theDOY = header.pFileCreationDate-... + datenum(... + str2double(... + datestr(header.pFileCreationDate,'yyyy')),1,1)+1; + end + + function theYear = get.fileCreationYear(header) + theYear = uint16(... + str2double(... + datestr(header.pFileCreationDate,'yyyy'))); + end + + function theVersion = get.version(header) + theVersion = [num2str(header.pVersionMajor) '.'... + num2str(header.pVersionMinor)]; + end + + function header = loadFrom(header, location) + % Load the LAS file public header from the given location. + + if (ischar(location)) + % Assume single argument is file name + [fid, msg]=fopen(location, 'r'); + + if (fid < 0) + error('LasPublicHeader:FileError',msg); + end + elseif (isnumeric(location)) + % Assume single argument is file ID + fid = location; + frewind(fid); + else + error('LasPublicHeader:InitError',... + 'Unknown argument initializer'); + end + + % Check the file signature to make sure we have a LAS file + if (~strcmp(sscanf(char(fread(fid,4,'uchar=>uchar')'),'%c'),... + header.fileSignature)) + error('LasPublicHeader:InvalidFile',... + 'File does not appear to be a valid LAS file.'); + end + + header.fileSourceID = fread(fid,1,'uint16=>uint16'); + header.globalEncoding = fread(fid,1,'uint16=>uint16'); + header.projectID_GUID1 = fread(fid,1,'uint32=>uint32'); + header.projectID_GUID2 = fread(fid,1,'uint16=>uint16'); + header.projectID_GUID3 = fread(fid,1,'uint16=>uint16'); + header.projectID_GUID4 = ... + sscanf(char(fread(fid,8,'uchar=>uchar')'),'%c'); + header.pVersionMajor = fread(fid,1,'uchar=>uchar'); + header.pVersionMinor = fread(fid,1,'uchar=>uchar'); + header.systemIdentifier= ... + sscanf(char(fread(fid,32,'uchar=>uchar')'),'%c'); + header.generatingSoftware = ... + sscanf(char(fread(fid,32,'uchar=>uchar')'),'%c'); + dayOfYear = fread(fid,1,'uint16=>double'); + year = fread(fid,1,'uint16=>double'); + header.pFileCreationDate = datenum(year,1,dayOfYear); + header.headerSize = fread(fid,1,'uint16=>uint16'); + header.offsetToPointData = fread(fid,1,'uint32=>uint32'); + header.nVLRs = fread(fid,1,'uint32=>uint32'); + header.pointDataFormatID = fread(fid,1,'uint8=>uint8'); + header.pointDataRecordLength = fread(fid,1,'uint16=>uint16'); + header.nPoints = fread(fid,1,'uint32=>uint32'); + header.nPointsByReturn = fread(fid,5,'uint32=>uint32'); + header.xScaleFactor = fread(fid,1,'double=>double'); + header.yScaleFactor = fread(fid,1,'double=>double'); + header.zScaleFactor = fread(fid,1,'double=>double'); + header.xOffset = fread(fid,1,'double=>double'); + header.yOffset = fread(fid,1,'double=>double'); + header.zOffset = fread(fid,1,'double=>double'); + header.maxX = fread(fid,1,'double=>double'); + header.minX = fread(fid,1,'double=>double'); + header.maxY = fread(fid,1,'double=>double'); + header.minY = fread(fid,1,'double=>double'); + header.maxZ = fread(fid,1,'double=>double'); + header.minZ = fread(fid,1,'double=>double'); + + if (header.headerSize > 227) + header.extraData = fread(fid,header.headerSize - 227,'uint8=>uint8'); + end + + if (ischar(location)) + % Close the file if it was opened here. + fclose(fid); + end + end + + function saveTo(header, location) + % Save the LAS file public header to the given location. + + if (ischar(location)) + % Assume single argument is file name + [fid, msg]=fopen(location, 'w'); + + if (fid < 0) + error('LasPublicHeader:FileError',msg); + end + elseif (isnumeric(location)) + % Assume single argument is file ID + fid = location; + else + error('LasPublicHeader:InitError',... + 'Unknown argument initializer'); + end + + fprintf(fid,'%4s',header.fileSignature); + fwrite(fid,header.fileSourceID,'uint16'); + fwrite(fid,header.globalEncoding,'uint16'); + fwrite(fid,header.projectID_GUID1,'uint32'); + fwrite(fid,header.projectID_GUID2,'uint16'); + fwrite(fid,header.projectID_GUID3,'uint16'); - % The Y Coordinate Scale Factor - % - % The scale factor fields contain a double floating point value - % that is uysed to scale the corresponding X, Y, and Z long values - % within the point records. The corresponding X, Y, and Z scale - % factor must be multiplied by the X, Y, or Z point record value to - % get the actual X, Y, or Z coordinate. For example, if the X, Y, - % and Z coordinates are intended to have two decimal point values, - % then each scale factor will contain the number 0.01. - yScaleFactor = double(0) + if (length(header.projectID_GUID4) > 8) + fwrite(fid,header.projectID_GUID4(1:8),'char'); + else + GUID4=[header.projectID_GUID4... + repmat(char(0),1,8-numel(header.projectID_GUID4))]; + fwrite(fid,GUID4,'char'); + end + + fprintf(fid,'%c',header.pVersionMajor); + fprintf(fid,'%c',header.pVersionMinor); - % The Z Coordinate Scale Factor - % - % The scale factor fields contain a double floating point value - % that is uysed to scale the corresponding X, Y, and Z long values - % within the point records. The corresponding X, Y, and Z scale - % factor must be multiplied by the X, Y, or Z point record value to - % get the actual X, Y, or Z coordinate. For example, if the X, Y, - % and Z coordinates are intended to have two decimal point values, - % then each scale factor will contain the number 0.01. - zScaleFactor = double(0) - - % The X Coordinate Offset - % - % The offset fields should be used to set the overall offset for - % the point records. In general these numbers will be zero, but for - % certain cases the resolution of the point data may not be large - % enough for a given projection system. However, it should always - % be assumed that these numbers are used. So to scale a given X - % from the point record, take the point record X multiplied by the - % X scale factor, and then add the X offset. - % - % X_coordinate = (X_record * X_scale) + X_offset - % Y_coordinate = (Y_record * Y_scale) + Y_offset - % Z_coordinate = (Z_record * Z_scale) + Z_offset - xOffset = double(0) - - % The Y Coordinate Offset - % - % The offset fields should be used to set the overall offset for - % the point records. In general these numbers will be zero, but for - % certain cases the resolution of the point data may not be large - % enough for a given projection system. However, it should always - % be assumed that these numbers are used. So to scale a given X - % from the point record, take the point record X multiplied by the - % X scale factor, and then add the X offset. - % - % X_coordinate = (X_record * X_scale) + X_offset - % Y_coordinate = (Y_record * Y_scale) + Y_offset - % Z_coordinate = (Z_record * Z_scale) + Z_offset - yOffset = double(0) - - % The Z Coordinate Offset - % - % The offset fields should be used to set the overall offset for - % the point records. In general these numbers will be zero, but for - % certain cases the resolution of the point data may not be large - % enough for a given projection system. However, it should always - % be assumed that these numbers are used. So to scale a given X - % from the point record, take the point record X multiplied by the - % X scale factor, and then add the X offset. - % - % X_coordinate = (X_record * X_scale) + X_offset - % Y_coordinate = (Y_record * Y_scale) + Y_offset - % Z_coordinate = (Z_record * Z_scale) + Z_offset - zOffset = double(0) - - % The Maximum X Coordinate - % - % The max and min data fields are the actual unscaled extents of - % the LAS point file data, specified in the coordinate system of - % the LAS data. - maxX = double(0) - - % The Minimum X Coordinate - % - % The max and min data fields are the actual unscaled extents of - % the LAS point file data, specified in the coordinate system of - % the LAS data. - minX = double(0) - - % The Maximum Y Coordinate - % - % The max and min data fields are the actual unscaled extents of - % the LAS point file data, specified in the coordinate system of - % the LAS data. - maxY = double(0) - - % The Minimum Y Coordinate - % - % The max and min data fields are the actual unscaled extents of - % the LAS point file data, specified in the coordinate system of - % the LAS data. - minY = double(0) - - % The Maximum Z Coordinate - % - % The max and min data fields are the actual unscaled extents of - % the LAS point file data, specified in the coordinate system of - % the LAS data. - maxZ = double(0) - - % The Minimum Z Coordinate - % - % The max and min data fields are the actual unscaled extents of - % the LAS point file data, specified in the coordinate system of - % the LAS data. - minZ = double(0) - - % The Extra Data - % - % The public header may be extended, though it is highly - % discouraged. If it is extended, we have no way of knowing how to - % interpret the extra data. We can simply pass it through, though. - extraData = zeros(0,'uint8'); - end - - properties (Access = private) - % The File Creation Date - % - % The date on which the file was created expressed as a serial date - % number. - pFileCreationDate - - % The Major Version Number - % - % The version number consists of a major and minor field. The major - % and minor fields combine to form the number that indicates the - % format number of the current specification itself. For example, - % specification number 1.2 (this version) would contain 1 in the - % major field and 2 in the minor field. - pVersionMajor = uint8(1) - - % The Minor Version Number - % - % The version number consists of a major and minor field. The major - % and minor fields combine to form the number that indicates the - % format number of the current specification itself. For example, - % specification number 1.2 (this version) would contain 1 in the - % major field and 2 in the minor field. - pVersionMinor = uint8(2) - end - - properties (Dependent) - % The File Creation Date - % - % The date on which the file was created expressed as a date - % string. - fileCreationDate - - % The File Creation Day of Year - % - % The day, expressed as an unsigned short, on which this file was - % created. Day is computed as the Greenwhich Mean Time (GMT) day. - % January 1 is considered day 1. - fileCreationDayOfYear - - % The File Creation Year - % - % The year, expressed as a four digit number, in which the file was - % created. - fileCreationYear - - % The Version Number - % - % The version number consists of a major and minor field. The major - % and minor fields combine to form the number that indicates the - % format number of the current specification itself. For example, - % specification number 1.2 (this version) would contain 1 in the - % major field and 2 in the minor field. - version - end - - methods - function header = LasPublicHeader(varargin) - % Constructor for LasFile public header block objects - % - % LasPublicHeader objects may be constructed either by using a - % no-argument (default) constructor or by passing a file - % descriptor or file name. - switch (nargin) - case 0 - header.pFileCreationDate = floor(datenum(1970,1,1)+... - java.lang.System.currentTimeMillis/... - (1e3*86400)); - case 1 - header = header.loadFrom(varargin{1}); - otherwise - error('LasPublicHeader:UnexpectedInputs',... - 'Unexpected number of inputs encountered.'); - end - end - - function theDate = get.fileCreationDate(header) - theDate = datestr(header.pFileCreationDate); - end - - function theDOY = get.fileCreationDayOfYear(header) - theDOY = header.pFileCreationDate-... - datenum(... - str2double(... - datestr(header.pFileCreationDate,'yyyy')),... - 1,1)+1; - end - - function theYear = get.fileCreationYear(header) - theYear = uint16(... - str2double(... - datestr(header.pFileCreationDate,'yyyy'))); - end - - function theVersion = get.version(header) - theVersion = [num2str(header.pVersionMajor) '.'... - num2str(header.pVersionMinor)]; - end - - function header = loadFrom(header, location) - % Load the LAS file public header from the given location. - - if (ischar(location)) - % Assume single argument is file name - [fid, msg]=fopen(location, 'r'); - - if (fid < 0) - error('LasPublicHeader:FileError',msg); - end - elseif (isnumeric(location)) - % Assume single argument is file ID - fid = location; - frewind(fid); - else - error('LasPublicHeader:InitError',... - 'Unknown argument initializer'); - end - - % Check the file signature to make sure we have a LAS file - if (~strcmp(sscanf(char(fread(fid,4,'uchar=>uchar')'),'%c'),... - header.fileSignature)) - error('LasPublicHeader:InvalidFile',... - 'File does not appear to be a valid LAS file.'); - end - - header.fileSourceID = fread(fid,1,'uint16=>uint16'); - header.globalEncoding = fread(fid,1,'uint16=>uint16'); - header.projectID_GUID1 = fread(fid,1,'uint32=>uint32'); - header.projectID_GUID2 = fread(fid,1,'uint16=>uint16'); - header.projectID_GUID3 = fread(fid,1,'uint16=>uint16'); - header.projectID_GUID4 = ... - sscanf(char(fread(fid,8,'uchar=>uchar')'),'%c'); - header.pVersionMajor = fread(fid,1,'uchar=>uchar'); - header.pVersionMinor = fread(fid,1,'uchar=>uchar'); - header.systemIdentifier= ... - sscanf(char(fread(fid,32,'uchar=>uchar')'),'%c'); - header.generatingSoftware = ... - sscanf(char(fread(fid,32,'uchar=>uchar')'),'%c'); - dayOfYear = fread(fid,1,'uint16=>double'); - year = fread(fid,1,'uint16=>double'); - header.pFileCreationDate = datenum(year,1,dayOfYear); - header.headerSize = fread(fid,1,'uint16=>uint16'); - header.offsetToPointData = fread(fid,1,'uint32=>uint32'); - header.nVLRs = fread(fid,1,'uint32=>uint32'); - header.pointDataFormatID = fread(fid,1,'uint8=>uint8'); - header.pointDataRecordLength = fread(fid,1,'uint16=>uint16'); - header.nPoints = fread(fid,1,'uint32=>uint32'); - header.nPointsByReturn = fread(fid,5,'uint32=>uint32'); - header.xScaleFactor = fread(fid,1,'double=>double'); - header.yScaleFactor = fread(fid,1,'double=>double'); - header.zScaleFactor = fread(fid,1,'double=>double'); - header.xOffset = fread(fid,1,'double=>double'); - header.yOffset = fread(fid,1,'double=>double'); - header.zOffset = fread(fid,1,'double=>double'); - header.maxX = fread(fid,1,'double=>double'); - header.minX = fread(fid,1,'double=>double'); - header.maxY = fread(fid,1,'double=>double'); - header.minY = fread(fid,1,'double=>double'); - header.maxZ = fread(fid,1,'double=>double'); - header.minZ = fread(fid,1,'double=>double'); - - if (header.headerSize > 227) - header.extraData = fread(fid,header.headerSize - 227,'uint8=>uint8'); - end - - if (ischar(location)) - % Close the file if it was opened here. - fclose(fid); - end - end - - function saveTo(header, fileid) - % Save the LAS file public header to the given location. - - if (ischar(fileid)) - % Assume single argument is file name - [fid, msg]=fopen(fileid, 'w'); - - if (fid < 0) - error('LasPublicHeader:FileError',msg); - end - elseif (isnumeric(fileid)) - % Assume single argument is file ID - fid = fileid; - else - error('LasPublicHeader:InitError',... - 'Unknown argument initializer'); - end - - fprintf(fid,'%4s',header.fileSignature); - fwrite(fid,header.fileSourceID,'uint16'); - fwrite(fid,header.globalEncoding,'uint16'); - fwrite(fid,header.projectID_GUID1,'uint32'); - fwrite(fid,header.projectID_GUID2,'uint16'); - fwrite(fid,header.projectID_GUID3,'uint16'); - if (length(header.projectID_GUID4) > 8) - fwrite(fid,header.projectID_GUID4(1:8),'char'); - else - GUID4=[header.projectID_GUID4... - repmat(char(0),1,8-numel(header.projectID_GUID4))]; - fwrite(fid,GUID4,'char'); - end - fprintf(fid,'%c',header.pVersionMajor); - fprintf(fid,'%c',header.pVersionMinor); - if (length(header.systemIdentifier) > 32) - fwrite(fid,header.systemIdentifier(1:32),'char'); - else - sysID=[header.systemIdentifier... - repmat(char(0),1,32-... - numel(header.systemIdentifier))]; - fwrite(fid,sysID,'char'); - end - if (length(header.generatingSoftware) > 32) - fwrite(fid,header.generatingSoftware(1:32),'char'); - else - genSoft=[header.generatingSoftware... - repmat(char(0),1,32-... - numel(header.generatingSoftware))]; - fwrite(fid,genSoft,'char'); - end - fwrite(fid,header.fileCreationDayOfYear,'uint16'); - fwrite(fid,header.fileCreationYear,'uint16'); - - fwrite(fid,header.headerSize,'uint16'); - fwrite(fid,header.offsetToPointData,'uint32'); - - fwrite(fid,header.nVLRs,'uint32'); - fwrite(fid,header.pointDataFormatID,'uint8'); - fwrite(fid,header.pointDataRecordLength,'uint16'); - fwrite(fid,header.nPoints,'uint32'); - fwrite(fid,header.nPointsByReturn,'uint32'); - fwrite(fid,header.xScaleFactor,'double'); - fwrite(fid,header.yScaleFactor,'double'); - fwrite(fid,header.zScaleFactor,'double'); - fwrite(fid,header.xOffset,'double'); - fwrite(fid,header.yOffset,'double'); - fwrite(fid,header.zOffset,'double'); - fwrite(fid,header.maxX,'double'); - fwrite(fid,header.minX,'double'); - fwrite(fid,header.maxY,'double'); - fwrite(fid,header.minY,'double'); - fwrite(fid,header.maxZ,'double'); - fwrite(fid,header.minZ,'double'); - fwrite(fid,header.extraData,'uint8'); - - if (ischar(fileid)) - % Close the file if it was opened here. - fclose(fid); - end - end - end + if (length(header.systemIdentifier) > 32) + fwrite(fid,header.systemIdentifier(1:32),'char'); + else + sysID=[header.systemIdentifier... + repmat(char(0),1,32-... + numel(header.systemIdentifier))]; + fwrite(fid,sysID,'char'); + end + + if (length(header.generatingSoftware) > 32) + fwrite(fid,header.generatingSoftware(1:32),'char'); + else + genSoft=[header.generatingSoftware... + repmat(char(0),1,32-... + numel(header.generatingSoftware))]; + fwrite(fid,genSoft,'char'); + end + + fwrite(fid,header.fileCreationDayOfYear,'uint16'); + fwrite(fid,header.fileCreationYear,'uint16'); + + fwrite(fid,header.headerSize,'uint16'); + fwrite(fid,header.offsetToPointData,'uint32'); + + fwrite(fid,header.nVLRs,'uint32'); + fwrite(fid,header.pointDataFormatID,'uint8'); + fwrite(fid,header.pointDataRecordLength,'uint16'); + fwrite(fid,header.nPoints,'uint32'); + fwrite(fid,header.nPointsByReturn,'uint32'); + fwrite(fid,header.xScaleFactor,'double'); + fwrite(fid,header.yScaleFactor,'double'); + fwrite(fid,header.zScaleFactor,'double'); + fwrite(fid,header.xOffset,'double'); + fwrite(fid,header.yOffset,'double'); + fwrite(fid,header.zOffset,'double'); + fwrite(fid,header.maxX,'double'); + fwrite(fid,header.minX,'double'); + fwrite(fid,header.maxY,'double'); + fwrite(fid,header.minY,'double'); + fwrite(fid,header.maxZ,'double'); + fwrite(fid,header.minZ,'double'); + fwrite(fid,header.extraData,'uint8'); + + if (ischar(location)) + % Close the file if it was opened here. + fclose(fid); + end + end + end end Index: Damkjer/Util/FileIO/LasFile/LasVlr.m =================================================================== --- Damkjer/Util/FileIO/LasFile/LasVlr.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/LasVlr.m (revision 6) @@ -1,179 +1,219 @@ -% LasVlr -% -% File: LasVlr.m -% -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% file variable length record. -% +% LasVlr ASPRS LAS format variable length record +% +% File: +% LasVlr.m +% +% Description: +% This MATLAB class represents an ASPRS LAS version 1.2 file variable length +% record. +% +% Limitations: +% None. +% +% Properties: +% userID - The User Identifier +% recordID - The Record Identifier +% recordLengthAfterHeader - The Record Length After Header +% description - The Description +% data - The Data +% +% Methods: +% [this] = LasVlr(varargin) - Constructor for LAS VLR objects. +% saveTo(location) - Save LAS VLR to a given location. +% +% Other m-files required: +% None. +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% +% See Also: +% LasFile +% + % Software History: -% 2012-AUG-29 K. Damkjer -% Initial Coding. +% 2012-AUG-29 K. Damkjer +% Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % classdef LasVlr - properties - % The User Identifier - % - % The User ID field is ASCII character data that identifies the - % user which created the variable length record. It is possible to - % have many Variable Length Records from different sources with - % different User IDs. If the character data is less than 16 - % characters, the remaining data must be null. The User ID must be - % registered with the LAS specification managing body. The - % management of these User IDs ensures that no two individuals - % accidentally use the same User ID. The specification will - % initially use two IDs: one for globally specified records - % (LASF_Spec), and another for projection types (LASF_Projection). - % Keys may be requested at - % http://www.asprs.org/lasform/keyform.html. - userID = '' - - % The Record Identifier - % - % The Record ID is dependent upon the User ID. There can be 0 to - % 65535 Record IDs for every User ID. The LAS specification manages - % its own Record IDs (User IDs owned by the specification), - % otherwise Record IDs will be managed by the owner of the given - % User ID. Thus each User ID is allowed to assign 0 to 65535 Record - % IDs in any manner they desire. Publicizing the meaning of a given - % Record ID is left to the owner of the given User ID. Unknown - % User ID/Record ID combinations should be ignored. - recordID = uint16(0) - - % The Record Length After Header - % - % The record length is the number of bytes for the record after the - % end of the standard part of the header. Thus the entire record - % length is 54 bytes (the header size iversion 1.2) plus the number - % of bytes in the variable length portion of the record. - recordLengthAfterHeader = uint16(0) - - % The Description - % - % Optional, null terminated description of the data. Any remaining - % characters not used must be null. - description = '' - - % The Data - % - % The raw VLR data as a byte array. - data = zeros(0,'uint8'); - end - - properties (Dependent) - end - - methods - function this = LasVlr(varargin) - % Constructor for LAS VLR Header objects + properties + % The User Identifier + % + % The User ID field is ASCII character data that identifies the + % user which created the variable length record. It is possible to + % have many Variable Length Records from different sources with + % different User IDs. If the character data is less than 16 + % characters, the remaining data must be null. The User ID must be + % registered with the LAS specification managing body. The + % management of these User IDs ensures that no two individuals + % accidentally use the same User ID. The specification will + % initially use two IDs: one for globally specified records + % (LASF_Spec), and another for projection types (LASF_Projection). + % Keys may be requested at + % http://www.asprs.org/lasform/keyform.html. + userID = '' + + % The Record Identifier + % + % The Record ID is dependent upon the User ID. There can be 0 to + % 65535 Record IDs for every User ID. The LAS specification manages + % its own Record IDs (User IDs owned by the specification), + % otherwise Record IDs will be managed by the owner of the given + % User ID. Thus each User ID is allowed to assign 0 to 65535 Record + % IDs in any manner they desire. Publicizing the meaning of a given + % Record ID is left to the owner of the given User ID. Unknown + % User ID/Record ID combinations should be ignored. + recordID = uint16(0) + + % The Record Length After Header + % + % The record length is the number of bytes for the record after the + % end of the standard part of the header. Thus the entire record + % length is 54 bytes (the header size iversion 1.2) plus the number + % of bytes in the variable length portion of the record. + recordLengthAfterHeader = uint16(0) + + % The Description + % + % Optional, null terminated description of the data. Any remaining + % characters not used must be null. + description = '' + + % The Data + % + % The raw VLR data as a byte array. + data = zeros(0,'uint8'); + end + + properties (Dependent) + end + + methods + function this = LasVlr(varargin) + % Constructor for LAS VLR objects + + switch (nargin) + case 0 + case 1 + if (isa(varargin{1}, 'LasVlr')) + % Copy constructor + rhs = varargin{1}; + fns = properties(rhs); + for i=1:length(fns) + this.(fns{i}) = rhs.(fns{i}); + end + else + if (ischar(varargin{1})) + % Assume single argument is file name + [fid, msg]=fopen(varargin{1}, 'rb'); + + if (fid < 0) + error('LasVlr:FileError',msg); + end + + % Check the file signature to make sure we have a LAS file + if (~strcmp(sscanf(char(fread(fid,4,'uchar=>uchar')'),... + '%c'),... + 'LASF')) + error('LasVlr:InvalidFile',... + 'File does not appear to be a valid LAS file.'); + end + elseif (isnumeric(varargin{1})) + % Assume single argument is file ID + fid = varargin{1}; + else + error('LasVlr:InitError',... + ['Unknown argument initializer: ' varargin{1}]); + end + + fseek(fid, 94, 'bof'); + headerSize = fread(fid,1,'uint16=>uint16'); + + fseek(fid, 100, 'bof'); + nVLRs = fread(fid,1,'uint16=>uint16'); + + this(nVLRs) = LasVlr; + + fseek(fid,headerSize,'bof'); + + for i=1:nVLRs + fread(fid,1,'uint16=>uint16'); + this(i).userID = ... + sscanf(char(fread(fid,16,'uchar=>uchar')'),'%c'); + this(i).recordID = ... + fread(fid,1,'uint16=>uint16'); + this(i).recordLengthAfterHeader = ... + fread(fid,1,'uint16=>uint16'); + this(i).description = ... + sscanf(char(fread(fid,32,'uchar=>uchar')'),'%c'); + this(i).data = fread(fid,... + this(i).recordLengthAfterHeader,'uint8=>uint8'); + end + + if (ischar(varargin{1})) + % Close the file if it was opened here. + fclose(fid); + end + end + otherwise + error('LasVlr:UnexpectedInputs',... + 'Unexpected number of inputs encountered.'); + end + end + + function saveTo(this, location) + % Save the LAS file public header to the given location. + + if (ischar(location)) + % Assume single argument is file name + [fid, msg]=fopen(location, 'w'); - switch (nargin) - case 0 - case 1 - if (ischar(varargin{1})) - % Assume single argument is file name - [fid, msg]=fopen(varargin{1}, 'rb'); - - if (fid < 0) - error('LasVlrHeader:FileError',msg); - end - - % Check the file signature to make sure we have a LAS file - if (~strcmp(sscanf(char(fread(fid,4,'uchar=>uchar')'),... - '%c'),... - 'LASF')) - error('LasVlrHeader:InvalidFile',... - 'File does not appear to be a valid LAS file.'); - end - elseif (isnumeric(varargin{1})) - % Assume single argument is file ID - fid = varargin{1}; - else - error('LasVlrHeader:InitError',... - 'Unknown argument initializer'); - end - - fseek(fid, 94, 'bof'); - headerSize = fread(fid,1,'uint16=>uint16'); - - fseek(fid, 100, 'bof'); - nVLRs = fread(fid,1,'uint16=>uint16'); - - this(nVLRs) = LasVlr; - - fseek(fid,headerSize,'bof'); - - for i=1:nVLRs - fread(fid,1,'uint16=>uint16'); - this(i).userID = ... - sscanf(char(fread(fid,16,'uchar=>uchar')'),'%c'); - this(i).recordID = ... - fread(fid,1,'uint16=>uint16'); - this(i).recordLengthAfterHeader = ... - fread(fid,1,'uint16=>uint16'); - this(i).description = ... - sscanf(char(fread(fid,32,'uchar=>uchar')'),'%c'); - this(i).data = fread(fid,... - this(i).recordLengthAfterHeader,'uint8=>uint8'); - end - - if (ischar(varargin{1})) - % Close the file if it was opened here. - fclose(fid); - end - otherwise - error('LasVlrHeader:UnexpectedInputs',... - 'Unexpected number of inputs encountered.'); + if (fid < 0) + error('LasPublicHeader:FileError',msg); end - end - - function saveTo(this, location) - % Save the LAS file public header to the given location. - - if (ischar(location)) - % Assume single argument is file name - [fid, msg]=fopen(location, 'w'); - - if (fid < 0) - error('LasPublicHeader:FileError',msg); - end - elseif (isnumeric(location)) - % Assume single argument is file ID - fid = location; - else - error('LasPublicHeader:InitError',... - 'Unknown argument initializer'); - end - - fwrite(fid,0,'uint16'); - - if (length(this.userID) > 16) - fwrite(fid,this.userID(1:16),'char'); - else - user=[this.userID... - repmat(char(0),1,16-... - numel(this.userID))]; - fwrite(fid,user,'char'); - end - - fwrite(fid,this.recordID,'uint16'); - fwrite(fid,this.recordLengthAfterHeader,'uint16'); - - if (length(this.description) > 32) - fwrite(fid,this.description(1:32),'char'); - else - desc=[this.description... - repmat(char(0),1,32-... - numel(this.description))]; - fwrite(fid,desc,'char'); - end - - fwrite(fid,this.data,'uint8'); - - if (ischar(location)) - % Close the file if it was opened here. - fclose(fid); - end - end - end + elseif (isnumeric(location)) + % Assume single argument is file ID + fid = location; + else + error('LasPublicHeader:InitError',... + 'Unknown argument initializer'); + end + + fwrite(fid,0,'uint16'); + + if (length(this.userID) > 16) + fwrite(fid,this.userID(1:16),'char'); + else + user=[this.userID... + repmat(char(0),1,16-... + numel(this.userID))]; + fwrite(fid,user,'char'); + end + + fwrite(fid,this.recordID,'uint16'); + fwrite(fid,this.recordLengthAfterHeader,'uint16'); + + if (length(this.description) > 32) + fwrite(fid,this.description(1:32),'char'); + else + desc=[this.description... + repmat(char(0),1,32-... + numel(this.description))]; + fwrite(fid,desc,'char'); + end + + fwrite(fid,this.data,'uint8'); + + if (ischar(location)) + % Close the file if it was opened here. + fclose(fid); + end + end + end end Index: Damkjer/Util/FileIO/LasFile/LasVlrOld.m =================================================================== --- Damkjer/Util/FileIO/LasFile/LasVlrOld.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/LasVlrOld.m (revision 6) @@ -1,18 +1,44 @@ -% LasfProjGeoAsciiParamsTagVLR +% LasVlrOld Obsolete Abstract Class for ASPRS LAS Format VLRs % -% File: LasfProjGeoAsciiParamsTagVLR.m +% File: +% LasVlrOld.m % -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% LASF Projection GeoAsciiParamsTag variable length record. +% Description: +% This class was intended to define and abstract interface for LAS VLRs. It +% is now obsolete. % +% Limitations: +% Obsolete Class. +% +% Properties: +% None. +% +% Methods: +% [vlr] = loadFrom(location) - Load LAS VLR from a gicen location. +% saveTo(location) - Save LAS VLR to a given location. +% +% Other m-files required: +% None. +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% +% See Also: +% LasFile + % Software History: % 2012-AUG-29 K. Damkjer % Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % classdef LasVlrOld - methods (Abstract = true) - saveTo(vlr,location) - vlr=loadFrom(vlr,location) - end + methods (Abstract = true) + saveTo(vlr,location) + vlr=loadFrom(vlr,location) + end end Index: Damkjer/Util/FileIO/LasFile/LasfProjGeoAsciiParamsTagVLR.m =================================================================== --- Damkjer/Util/FileIO/LasFile/LasfProjGeoAsciiParamsTagVLR.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/LasfProjGeoAsciiParamsTagVLR.m (revision 6) @@ -1,32 +1,81 @@ -% LasfProjGeoAsciiParamsTagVLR +% LasfProjGeoAsciiParamsTagVLR LASF_Projection: GeoTIFF Tag - ASCII Parameters % -% File: LasfProjGeoAsciiParamsTagVLR.m +% File: +% LasfProjGeoAsciiParamsTagVLR.m % -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% LASF Projection GeoAsciiParamsTag variable length record. +% Description: +% This MATLAB class represents a VLR that contains an array of ASCII data. It +% contains many strings separated by null terminator characters which are +% referenced by position from data in the GeoKeyDirectoryTag record. % +% Limitations: +% None. +% +% Properties: +% params - The ASCII Parameters +% userID - The User Identifier +% recordID - The Record Identifier +% recordLengthAfterHeader - The Record Length After Header +% description - The Description +% data - The Data +% +% Methods: +% [this] = LasfProjGeoAsciiParamsTagVLR(varargin) - Constructor. +% saveTo(location) - Save LAS VLR to a given +% location. +% +% Other m-files required: +% LasVlr.m +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% http://www.remotesensing.org/geotiff/spec/geotiffhome.html +% +% See Also: +% LasFile +% LasVlr +% LasfProjGeoKeyDirectoryTagVLR +% + % Software History: % 2012-AUG-29 K. Damkjer % Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % classdef (Sealed = true) LasfProjGeoAsciiParamsTagVLR < LasVlr - properties - params - end - - methods - function vlr = LasfProjGeoAsciiParamsTagVLR(varargin) - if (nargin == 1) - if (~strcmp(class(varargin{1}),'uint8')) - error('LasfProjGeoAsciiParamsTagVLR:InitError',... - 'Unknown argument initializer'); - end - - % The GeoKeyDirectoryTag data is a stream of unsigned - % shorts - vlr.params = char(varargin{1}); + properties + params + end + + methods + function vlr = LasfProjGeoAsciiParamsTagVLR(varargin) + if (nargin == 0) + superargs = {}; + elseif (nargin == 1) + if (~isa(varargin{1},'LasVlr')) + error('LasfProjGeoDoubleParamsTagVLR:InitError',... + 'Unknown argument initializer.'); end - end - end + + vlrHeader = varargin{1}; + + superargs{1} = vlrHeader; + else + error('LasfProjGeoDoubleParamsTagVLR:InitError',... + 'Unexpected number of inputs.'); + end + + vlr = vlr@LasVlr(superargs{:}); + + if (nargin == 1) + % The GeoAsciiParamsTag data is a stream of null-terminated strings + vlr.params = char(vlrHeader.data); + end + end + end end Index: Damkjer/Util/FileIO/LasFile/LasfProjGeoDoubleParamsTagVLR.m =================================================================== --- Damkjer/Util/FileIO/LasFile/LasfProjGeoDoubleParamsTagVLR.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/LasfProjGeoDoubleParamsTagVLR.m (revision 6) @@ -1,32 +1,80 @@ -% LasfProjGeoAsciiParamsTagVLR +% LasfProjGeoDoubleParamsTagVLR LASF_Projection: GeoTIFF Tag - Double Parameters % -% File: LasfProjGeoAsciiParamsTagVLR.m +% File: +% LasfProjGeoDoubleParamsTagVLR.m % -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% LASF Projection GeoAsciiParamsTag variable length record. +% Description: +% This MATLAB class represents a VLR that contains an array of double values +% referenced by tag sets in the GeoKeyDirectoryTag record. % +% Limitations: +% None. +% +% Properties: +% params - The Double Parameters +% userID - The User Identifier +% recordID - The Record Identifier +% recordLengthAfterHeader - The Record Length After Header +% description - The Description +% data - The Data +% +% Methods: +% [this] = LasfProjGeoDoubleParamsTagVLR(varargin) - Constructor. +% saveTo(location) - Save LAS VLR to a given +% location. +% +% Other m-files required: +% LasVlr.m +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% http://www.remotesensing.org/geotiff/spec/geotiffhome.html +% +% See Also: +% LasFile +% LasVlr +% LasfProjGeoKeyDirectoryTagVLR +% + % Software History: % 2012-AUG-29 K. Damkjer % Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % classdef (Sealed = true) LasfProjGeoDoubleParamsTagVLR < LasVlr - properties - params - end - - methods - function vlr = LasfProjGeoDoubleParamsTagVLR(varargin) - if (nargin == 1) - if (~strcmp(class(varargin{1}),'uint8')) - error('LasfProjGeoDoubleParamsTagVLR:InitError',... - 'Unknown argument initializer'); - end - - % The GeoKeyDirectoryTag data is a stream of unsigned - % shorts - vlr.params = typecast(varargin{1},'double'); + properties + params + end + + methods + function vlr = LasfProjGeoDoubleParamsTagVLR(varargin) + if (nargin == 0) + superargs = {}; + elseif (nargin == 1) + if (~isa(varargin{1},'LasVlr')) + error('LasfProjGeoDoubleParamsTagVLR:InitError',... + 'Unknown argument initializer.'); end - end - end + + vlrHeader = varargin{1}; + + superargs{1} = vlrHeader; + else + error('LasfProjGeoDoubleParamsTagVLR:InitError',... + 'Unexpected number of inputs.'); + end + + vlr = vlr@LasVlr(superargs{:}); + + if (nargin == 1) + % The GeoDoubleParamsTag data is a stream of doubles + vlr.params = typecast(vlrHeader.data,'double'); + end + end + end end Index: Damkjer/Util/FileIO/LasFile/LasfProjGeoKeyDirectoryTagVLR.m =================================================================== --- Damkjer/Util/FileIO/LasFile/LasfProjGeoKeyDirectoryTagVLR.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/LasfProjGeoKeyDirectoryTagVLR.m (revision 6) @@ -1,49 +1,117 @@ -% LasfProjGeoAsciiParamsTagVLR +% LasfProjGeoKeyDirectoryTagVLR LASF_Projection: GeoTIFF Tag - Key Directory % -% File: LasfProjGeoAsciiParamsTagVLR.m +% File: +% LasfProjGeoKeyDirectoryTagVLR.m % -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% LASF Projection GeoAsciiParamsTag variable length record. +% Description: +% This MATLAB class represents a LAS File Projection VLR containing the key +% values that define the coordinate system. A complete description can be +% found in the GeoTIFF format specification. % +% Limitations: +% None. +% +% Properties: +% keyDirectoryVersion - The Key Directory Version (1) +% keyRevision - The Key Revision (1) +% minorRevision - The Minor Revision (0) +% numberOfKeys - The Number of GeoTIFF Keys in this Directory +% keys - The GeoTIFF Keys +% userID - The User Identifier +% recordID - The Record Identifier +% recordLengthAfterHeader - The Record Length After Header +% description - The Description +% data - The Data +% +% Methods: +% [this] = LasfProjGeoKeyDirectoryTagVLR(varargin) - Constructor. +% saveTo(location) - Save LAS VLR to a given +% location. +% +% Other m-files required: +% LasVlr.m +% geoKeyFactory.m +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% http://www.remotesensing.org/geotiff/spec/geotiffhome.html +% +% See Also: +% LasFile +% LasVlr +% geoKeyFactory +% + % Software History: % 2012-AUG-29 K. Damkjer % Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % classdef (Sealed = true) LasfProjGeoKeyDirectoryTagVLR < LasVlr - properties - keyDirectoryVersion = 1 - keyRevision = 1 - minorRevision = 0 - numberOfKeys = 0 - end - - methods - function vlr = LasfProjGeoKeyDirectoryTagVLR(varargin) - if (nargin == 1) - if (~strcmp(class(varargin{1}),'uint8')) - error('LasfProjGeoKeyDirectoryTagVLR:InitError',... - 'Unknown argument initializer'); - end + properties + keyDirectoryVersion = 1 + keyRevision = 1 + minorRevision = 0 + numberOfKeys = 0 + keys = {} + end + + methods + function vlr = LasfProjGeoKeyDirectoryTagVLR(varargin) + if (nargin == 0) + superargs = {}; + elseif (nargin == 1) + if (~isa(varargin{1},'LasVlr')) + error('LasfProjGeoKeyDirectoryTagVLR:InitError',... + 'Unknown argument initializer.'); + end + + vlrHeader = varargin{1}; + + superargs{1} = vlrHeader; + else + error('LasfProjGeoKeyDirectoryTagVLR:InitError',... + 'Unexpected number of inputs.'); + end + + vlr = vlr@LasVlr(superargs{:}); + + if (nargin == 1) + + if (~(strcmp(vlrHeader.userID, 'LASF_Projection') &&... + vlrHeader.recordID == 34735)) + error('LasfProjGeoKeyDirectoryTagVLR:InitError',... + 'Invalid argument.'); + end + + % The GeoKeyDirectoryTag data is a stream of unsigned shorts + data = typecast(vlrHeader.data,'uint16'); + + if (data(1) ~= 1 && data(2) ~= 1 && data(3) ~= 0) + error('LasfProjGeoKeyDirectoryTagVLR:InitError',... + 'Invalid Format.'); + end + + vlr.numberOfKeys = data(4); + + keyData = reshape(data(5:end),4,vlr.numberOfKeys); - % The GeoKeyDirectoryTag data is a stream of unsigned - % shorts - data = typecast(varargin{1},'uint16'); - - if (data(1) ~= 1 && data(2) ~= 1 && data(3) ~= 0) - error('LasfProjGeoKeyDirectoryTagVLR:InitError',... - 'Invalid Format.'); - end - - vlr.numberOfKeys = data(4); - - keys = reshape(data(5:end),4,vlr.numberOfKeys); - - for key=1:vlr.numberOfKeys - tempKey=keys(:,key); - geoKeyFactory(tempKey(1),tempKey(2),tempKey(3),tempKey(4)) - end + vlr.keys = cell(vlr.numberOfKeys,1); + + for key=1:vlr.numberOfKeys + tempKey=keyData(:,key); + + if (~(tempKey(1) == 0 && tempKey(2) == 0)) + vlr.keys{key}=... + geoKeyFactory(tempKey(1),tempKey(2),tempKey(3),tempKey(4)); + end end - end - end + end + end + end end Index: Damkjer/Util/FileIO/LasFile/LasfSpecClassificationVLR.m =================================================================== --- Damkjer/Util/FileIO/LasFile/LasfSpecClassificationVLR.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/LasfSpecClassificationVLR.m (revision 6) @@ -1,51 +1,99 @@ +% LasfSpecClassificationVLR LASF_Spec: Classification Dictionary +% +% File: +% LasfSpecClassificationVLR.m +% +% Description: +% This MATLAB class represents an ASPRS LAS version 1.2 file variable length +% record. +% +% Limitations: +% None. +% +% Properties: +% classes - The Classification Dictionary +% userID - The User Identifier +% recordID - The Record Identifier +% recordLengthAfterHeader - The Record Length After Header +% description - The Description +% data - The Data +% +% Methods: +% [this] = LasfSpecClassificationVLR(varargin) - Constructor. +% saveTo(location) - Save LAS VLR to a given +% location. +% +% Other m-files required: +% LasVlr.m +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% +% See Also: +% LasFile +% LasVlr +% + +% Software History: +% 2012-AUG-29 K. Damkjer +% Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. +% + classdef (Sealed = true) LasfSpecClassificationVLR < LasVlr - properties (SetAccess = private) - classes = containers.Map(... - {0,... - 1,... - 2,... - 3,... - 4,... - 5,... - 6,... - 7,... - 8,... - 9,... - 12},... - {'Not Classified',... - 'Unclassified',... - 'Ground',... - 'Low Veg',... - 'Medium Veg',... - 'High Veg',... - 'Building',... - 'Noise',... - 'Mass Point',... - 'Water',... - 'Overlap Point'}); - end - - methods - function vlr = LasfSpecClassificationVLR(varargin) - if (nargin == 1) - if (~strcmp(class(varargin{1}),'uint8')) - error('LasfSpecClassificationVLR:InitError',... - 'Unknown argument initializer'); - end - - % The LasfSpecClassificationVLR data is a stream of uint8 - % char[15]s - data = reshape(varargin{1},16,length(varargin{1})/16); - - keySet=num2cell(data(1,:)); - valueSet=cell(1,size(data,2)); - - for i=1:size(data,2) - valueSet{i}=sscanf(char(data(2:end,i)),'%c'); - end - - vlr.classes = containers.Map(keySet,valueSet,'UniformValues',true); + properties (SetAccess = private) + classes = containers.Map(... + { 0,... + 1,... + 2,... + 3,... + 4,... + 5,... + 6,... + 7,... + 8,... + 9,... + 12},... + {'Not Classified',... + 'Unclassified',... + 'Ground',... + 'Low Veg',... + 'Medium Veg',... + 'High Veg',... + 'Building',... + 'Noise',... + 'Mass Point',... + 'Water',... + 'Overlap Point'}); + end + + methods + function vlr = LasfSpecClassificationVLR(varargin) + vlr = vlr@LasVlr(varargin{1}); + + if (nargin == 1) + if (~strcmp(class(varargin{1}),'uint8')) + error('LasfSpecClassificationVLR:InitError',... + 'Unknown argument initializer'); end - end - end + + % The LasfSpecClassificationVLR data is a stream of uint8 + % char[15]s + data = reshape(varargin{1},16,length(varargin{1})/16); + + keySet=num2cell(data(1,:)); + valueSet=cell(1,size(data,2)); + + for i=1:size(data,2) + valueSet{i}=sscanf(char(data(2:end,i)),'%c'); + end + + vlr.classes = containers.Map(keySet,valueSet,'UniformValues',true); + end + end + end end Index: Damkjer/Util/FileIO/LasFile/lasVlrFactory.m =================================================================== --- Damkjer/Util/FileIO/LasFile/lasVlrFactory.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/lasVlrFactory.m (revision 6) @@ -1,24 +1,59 @@ -% LasfProjGeoAsciiParamsTagVLR +% lasVlrFactory Abstract factory for LAS VLR objects % -% File: LasfProjGeoAsciiParamsTagVLR.m +% File: +% lasVlrFactory.m % -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% LASF Projection GeoAsciiParamsTag variable length record. +% Description: +% This MATLAB function provides an interface similar to an abstract factory +% design pattern for constructing known LAS VLR objects. +% +% Limitations: +% None. +% +% Synopsis: +% [vlr] = lasVlrFactory(vlrHeader) +% +% Inputs: +% vlrHeader - A LAS VLR object containing populated header data +% +% Outputs: +% vlr - A specialized LAS VLR object +% +% Other m-files required: +% LasVlr.m +% lasfProjVlrFactory.m +% lasfSpecVlrFactory.m +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% +% See Also: +% LasFile +% LasVlr +% lasfProjVlrFactory +% lasfSpecVlrFactory +% + +% Copyright (C) 2012 Kristian L. Damkjer. % % Software History: % 2012-AUG-29 K. Damkjer % Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % function [vlr] = lasVlrFactory(vlrHeader) switch vlrHeader.userID - case 'LASF_Projection' - vlr = lasfProjVlrFactory(vlrHeader.recordID, vlrHeader.data); - case 'LASF_Spec' - vlr = lasfSpecVlrFactory(vlrHeader.recordID, vlrHeader.data); - otherwise - disp(vlrHeader.userID); - disp(vlrHeader.description); -% vlr = LasUnknownVlr(vlrHeader); + case 'LASF_Projection' + vlr = lasfProjVlrFactory(vlrHeader); + case 'LASF_Spec' + vlr = lasfSpecVlrFactory(vlrHeader); + otherwise + disp(vlrHeader.userID); + disp(vlrHeader.description); end end Index: Damkjer/Util/FileIO/LasFile/lasfProjVlrFactory.m =================================================================== --- Damkjer/Util/FileIO/LasFile/lasfProjVlrFactory.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/lasfProjVlrFactory.m (revision 6) @@ -1,24 +1,64 @@ -% LasfProjGeoAsciiParamsTagVLR +% lasfProjVlrFactory Concrete factory for LAS File Projection VLR objects % -% File: LasfProjGeoAsciiParamsTagVLR.m +% File: +% lasfProjVlrFactory.m % -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% LASF Projection GeoAsciiParamsTag variable length record. +% Description: +% This MATLAB function provides an interface similar to a concrete factory +% in the abstract factory design pattern for constructing defined LAS File +% Projection VLR objects. % +% Limitations: +% None. +% +% Synopsis: +% [vlr] = lasfProjVlrFactory(vlrHeader) +% +% Inputs: +% vlrHeader - A LAS VLR object containing populated header data +% +% Outputs: +% vlr - A specialized LAS File VLR Projection object +% +% Other m-files required: +% LasVlr.m +% LasfProjGeoKeyDirectoryTagVLR.m +% LasfProjGeoDoubleParamsTagVLR.m +% LasfProjGeoAsciiParamsTagVLR.m +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% +% See Also: +% LasFile +% LasVlr +% lasVlrFactory +% LasfProjGeoKeyDirectoryTagVLR +% LasfProjGeoDoubleParamsTagVLR +% LasfProjGeoAsciiParamsTagVLR +% + % Software History: % 2012-AUG-29 K. Damkjer % Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % -function [vlr] = lasfProjVlrFactory(recordID, data) -switch recordID - case 34735 - vlr = LasfProjGeoKeyDirectoryTagVLR(data); - case 34736 - vlr = LasfProjGeoDoubleParamsTagVLR(data); - case 34737 - vlr = LasfProjGeoAsciiParamsTagVLR(data); - otherwise - vlr = LasUnknownVlr(data); +function [vlr] = lasfProjVlrFactory(vlrHeader) +switch vlrHeader.recordID + case 34735 + vlr = LasfProjGeoKeyDirectoryTagVLR(vlrHeader); + case 34736 + vlr = LasfProjGeoDoubleParamsTagVLR(vlrHeader); + case 34737 + vlr = LasfProjGeoAsciiParamsTagVLR(vlrHeader); + otherwise + disp(vlrHeader.userID); + disp(vlrHeader.recordID); + disp(vlrHeader.description); end end Index: Damkjer/Util/FileIO/LasFile/lasfSpecVlrFactory.m =================================================================== --- Damkjer/Util/FileIO/LasFile/lasfSpecVlrFactory.m (revision 0) +++ Damkjer/Util/FileIO/LasFile/lasfSpecVlrFactory.m (revision 6) @@ -1,24 +1,60 @@ -% LasfProjGeoAsciiParamsTagVLR +% lasfSpecVlrFactory Concrete factory for LAS File Specification VLR objects % -% File: LasfProjGeoAsciiParamsTagVLR.m +% File: +% lasfSpecVlrFactory.m % -% Description: This MATLAB class represents an ASPRS LAS version 1.2 -% LASF Projection GeoAsciiParamsTag variable length record. +% Description: +% This MATLAB function provides an interface similar to a concrete factory +% in the abstract factory design pattern for constructing defined LAS File +% Specification VLR objects. % +% Limitations: +% LasfSpecHistogramVLR and LasfSpecTextAreaVLR are currently undefined. +% +% Synopsis: +% [vlr] = lasfSpecVlrFactory(vlrHeader) +% +% Inputs: +% vlrHeader - A LAS VLR object containing populated header data +% +% Outputs: +% vlr - A specialized LAS File VLR Specification object +% +% Other m-files required: +% LasVlr.m +% LasfSpecClassification.m +% +% MAT-files required: +% None. +% +% References: +% http://asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf +% +% See Also: +% LasFile +% LasVlr +% lasVlrFactory +% LasfSpecClassificationVLR +% + % Software History: % 2012-AUG-29 K. Damkjer % Initial Coding. +% 2013-JUN-17 K. Damkjer +% Additional Commenting. % -function [vlr] = lasfSpecVlrFactory(recordID, data) -switch recordID +function [vlr] = lasfSpecVlrFactory(vlrHeader) +switch vlrHeader.recordID case 0 - vlr = LasfSpecClassificationVLR(data); + vlr = LasfSpecClassificationVLR(vlrHeader.data); case 2 - vlr = LasfSpecHistogramVLR(data); + vlr = LasfSpecHistogramVLR(vlrHeader.data); case 3 - vlr = LasfSpecTextAreaVLR(data); + vlr = LasfSpecTextAreaVLR(vlrHeader.data); otherwise - vlr = LasUnknownVlr(data); + disp(vlrHeader.userID); + disp(vlrHeader.recordID); + disp(vlrHeader.description); end end Index: Damkjer/Util/FileIO/PlyFile/readPly.m =================================================================== --- Damkjer/Util/FileIO/PlyFile/readPly.m (revision 0) +++ Damkjer/Util/FileIO/PlyFile/readPly.m (revision 6) @@ -287,7 +287,9 @@ propnames=fieldnames(element); - msg=['Loading element ' ielem ' of ' length(elemnames)]; + msg=['Loading element ' num2str(ielem) ' of '... + num2str(length(elemnames)) ' ']; tstart=tic; - h = timebar(1, elements, msg, tstart); + h = timebar(1, length(elemnames), msg, tstart); + tic; if (element.hasListProp) Index: Damkjer/Util/MATLAB/ClassHandle.cpp =================================================================== --- Damkjer/Util/MATLAB/ClassHandle.cpp (revision 0) +++ Damkjer/Util/MATLAB/ClassHandle.cpp (revision 6) @@ -1,2 +1,27 @@ +//************************************************************************* +// FILE: ClassHandle.cpp +// +// Copyright (C) 2012 Kristian Damkjer. +// +// DESCRIPTION: This class is intended to be used to wrap C++ objects for +// transport across mex function calls in MATLAB. +// +// LIMITATIONS: To preserve the lifetime guarantees of the pointer, default +// construction, copy construction and assignment are +// explicitly disallowed. +// +// The class handle owns the pointer it contains. Thus, when +// the class handle is destroyed, the pointer is deleted. +// +// Class handles are forced to be allocated on the heap by +// using the named constructor idiom. Class handles should be +// destroyed using the destroyHandleTo template function. +// +// SOFTWARE HISTORY: +//> 2012-OCT-08 K. Damkjer +// Initial Coding. +//< +//************************************************************************* + #include "ClassHandle.h" Index: Damkjer/Util/MATLAB/ClassHandle.h =================================================================== --- Damkjer/Util/MATLAB/ClassHandle.h (revision 0) +++ Damkjer/Util/MATLAB/ClassHandle.h (revision 6) @@ -2,5 +2,5 @@ // FILE: ClassHandle.h // -// Copyright (C) 2012 Kristian Damkjer. +// Copyright (C) 2012 Kristian L. Damkjer. // // DESCRIPTION: This class is intended to be used to wrap C++ objects for @@ -96,5 +96,5 @@ { mexErrMsgIdAndTxt("Damkjer:matAsPtr:invalidHandle", - "Input must be real-valued index-class scalar."); + "Input must be real-valued index-class scalar."); } @@ -122,5 +122,5 @@ { mexErrMsgIdAndTxt("Damkjer:matAsPtr:invalidHandle", - "Input must be real-valued index-class scalar."); + "Input must be real-valued index-class scalar."); } Index: Damkjer/Util/Math/fastcov.cpp =================================================================== --- Damkjer/Util/Math/fastcov.cpp (revision 3) +++ Damkjer/Util/Math/fastcov.cpp (revision 6) @@ -131,14 +131,14 @@ { // Center to mean. -// cov_vals[cellp][n2 + Ns[cellp] * n1] += -// w2 -// * (vals[cellp][mc + Ms[cellp] * n1]-mean[n1]) -// * (vals[cellp][mc + Ms[cellp] * n2]-mean[n2]); + cov_vals[cellp][n2 + Ns[cellp] * n1] += + w2 + * (vals[cellp][mc + Ms[cellp] * n1]-mean[n1]) + * (vals[cellp][mc + Ms[cellp] * n2]-mean[n2]); // Center to first point. - cov_vals[cellp][n2 + Ns[cellp] * n1] += - w2 - * (vals[cellp][mc + Ms[cellp] * n1]-vals[cellp][Ms[cellp] * n1]) - * (vals[cellp][mc + Ms[cellp] * n2]-vals[cellp][Ms[cellp] * n2]); +// cov_vals[cellp][n2 + Ns[cellp] * n1] += +// w2 +// * (vals[cellp][mc + Ms[cellp] * n1]-vals[cellp][Ms[cellp] * n1]) +// * (vals[cellp][mc + Ms[cellp] * n2]-vals[cellp][Ms[cellp] * n2]); // End comment. } Index: Damkjer/Util/Math/makeMath.m =================================================================== --- Damkjer/Util/Math/makeMath.m (revision 3) +++ Damkjer/Util/Math/makeMath.m (revision 6) @@ -25,48 +25,51 @@ flags = {'verbose', 'debug', 'warnings', 'parallel'}; -options=[]; +options = ''; while ~isempty(varargin) - flag = varargin{1}; - % If the property has been supplied in a shortened form, lengthen it - iFlag = find(strncmpi(flag, flags, length(flag))); - - if isempty(iFlag) - options=[options ' ' flag]; - error('Damkjer:makeVpsTree:InvalidFlag', ['Invalid Flag: ' flag]); - elseif length(iFlag) > 1 - error('Damkjer:makeColorize:AmbiguousFlag', ... - ['Supplied shortened flag is ambiguous: ' flag]); - else - flag = flags{iFlag}; - - switch flag - case 'verbose' - command = [command ' -v']; - case 'debug' - command = [command ' -g']; - case 'parallel' - if (ispc) - command = [command ' COMPFLAGS="$COMPFLAGS /openmp"']; - else - command = [command ' CFLAGS="\$CFLAGS -fopenmp"' ... - ' CXXFLAGS="\$CXXFLAGS -fopenmp"' ... - ' LDFLAGS="\$LDFLAGS -fopenmp"']; - end - case 'warnings' - if (ispc) - command = [command ' -v COMPFLAGS="$COMPFLAGS -W3"']; - else - command = [command ' -v' ... - ' CFLAGS="\$CFLAGS -Wall -Wextra"' ... - ' CXXFLAGS="\$CXXFLAGS -Wall -Wextra"']; - end - otherwise - error('Damkjer:makeVpsTree:NotImplementedFlag', ... - ['Flag recognized, but not supported (oops): ' flag]); - end % switch property - end - - varargin(1) = []; + flag = varargin{1}; + % If the property has been supplied in a shortened form, lengthen it + iFlag = find(strncmpi(flag, flags, length(flag))); + + if isempty(iFlag) + options = strcat(options, ' ', flag); + + % Pass through any arguments we don't recognize, alternatively escalate + % to an error. +% error('Damkjer:makeVpsTree:InvalidFlag', ['Invalid Flag: ' flag]); + elseif length(iFlag) > 1 + error('Damkjer:makeMath:AmbiguousFlag', ... + ['Supplied shortened flag is ambiguous: ' flag]); + else + flag = flags{iFlag}; + + switch flag + case 'verbose' + command = strcat(command, ' -v'); + case 'debug' + command = strcat(command, ' -g'); + case 'parallel' + if (ispc) + command = strcat(command, ' COMPFLAGS="$COMPFLAGS /openmp"'); + else + command = strcat(command, ' CFLAGS="\$CFLAGS -fopenmp"',... + ' CXXFLAGS="\$CXXFLAGS -fopenmp"',... + ' LDFLAGS="\$LDFLAGS -fopenmp"'); + end + case 'warnings' + if (ispc) + command = strcat(command, ' -v COMPFLAGS="$COMPFLAGS -W3"'); + else + command = strcat(command, ' -v',... + ' CFLAGS="\$CFLAGS -Wall -Wextra"',... + ' CXXFLAGS="\$CXXFLAGS -Wall -Wextra"'); + end + otherwise + error('Damkjer:makeVpsTree:NotImplementedFlag', ... + ['Flag recognized, but not supported (oops): ' flag]); + end % switch property + end + + varargin(1) = []; end % while Index: Damkjer/Util/Math/quadricfit.m =================================================================== --- Damkjer/Util/Math/quadricfit.m (revision 6) +++ Damkjer/Util/Math/quadricfit.m (revision 6) @@ -0,0 +1,37 @@ +function [ Q, center ] = quadricfit( pnts ) +%ELLIPSOIDFIT Summary of this function goes here +% Detailed explanation goes here + +size(pnts) + +if (size(pnts,1)==3) + pnts=[pnts;ones(1,size(pnts,2))]; +end + +if (size(pnts,1)~=4) + error 'Points must be 3D' +end + +A=[(pnts(1,:).^2)', ... %x^2 + (pnts(1,:).*pnts(2,:))', ... %xy + (pnts(1,:).*pnts(3,:))', ... %xz + (pnts(2,:).^2)', ... %y^2 + (pnts(2,:).*pnts(3,:))', ... %yz + (pnts(3,:).^2)', ... %z^2 + pnts(1,:)', ... %x + pnts(2,:)', ... %y + pnts(3,:)', ... %z + pnts(4,:)']; %w + +QQ=null(A); + +Q=[QQ(1), QQ(2)/2,QQ(3)/2,QQ(7)/2; + QQ(2)/2,QQ(4), QQ(5)/2,QQ(8)/2; + QQ(3)/2,QQ(5)/2,QQ(6), QQ(9)/2; + QQ(7)/2,QQ(8)/2,QQ(9)/2,QQ(10)]; + +Q=Q/Q(4,4); + +center=Q(1:3,1:3)\Q(1:3,4); +end + Index: Damkjer/Util/Math/timeFormat.m =================================================================== --- Damkjer/Util/Math/timeFormat.m (revision 6) +++ Damkjer/Util/Math/timeFormat.m (revision 6) @@ -0,0 +1,60 @@ +% timeFormat Format seconds as hh:mm:ss.000 +% +% File: +% timeFormat.m +% +% Description: +% This MATLAB function converts a number representing seconds to a string +% format as hh:mm:ss.000. +% +% Limitations: +% None. +% +% Synopsis: +% [msg] = timeFormat(seconds) +% +% Inputs: +% seconds - A numeric representing time in seconds. +% +% Outputs: +% msg - A string representation of seconds in hh:mm:ss.000 format. +% +% Toolbox requirements: +% None. +% +% Script requirements: +% None. +% +% Data requirements: +% None. +% +% References: +% None. +% +% See Also: +% None. +% + +% Copyright (C) 2013 Kristian L. Damkjer. +% +% Software History: +% 2013-JUN-20 K. Damkjer +% Initial Coding. +% + +function [msg] = timeFormat(seconds) + +h = floor(seconds/3600); +m = floor((seconds - h * 3600)/60); +s = mod(seconds, 60); + +if (h > 0) + msg = [num2str(h,'%02d') ':' num2str(m,'%02d') ':' num2str(s,'%06.3f')]; +elseif (m > 0) + msg = [num2str(m,'%02d') ':' num2str(s,'%06.3f')]; +else + msg = ['0:' num2str(s,'%06.3f')]; +end + +end + Index: Damkjer/Util/SpatialIndexing/VpTree/makeVpTree.m =================================================================== --- Damkjer/Util/SpatialIndexing/VpTree/makeVpTree.m (revision 0) +++ Damkjer/Util/SpatialIndexing/VpTree/makeVpTree.m (revision 6) @@ -25,52 +25,56 @@ flags = {'verbose', 'debug', 'warnings', 'parallel'}; -options=[]; +options = ''; while ~isempty(varargin) - flag = varargin{1}; - % If the property has been supplied in a shortened form, lengthen it - iFlag = find(strncmpi(flag, flags, length(flag))); - - if isempty(iFlag) - options=[options ' ' flag]; - error('Damkjer:makeVpsTree:InvalidFlag', ['Invalid Flag: ' flag]); - elseif length(iFlag) > 1 - error('Damkjer:makeVpsTree:AmbiguousFlag', ... - ['Supplied shortened flag is ambiguous: ' flag]); - else - flag = flags{iFlag}; - - switch flag - case 'verbose' - command = [command ' -v']; - case 'debug' - command = [command ' -g']; - case 'parallel' - if (ispc) - command = [command ' COMPFLAGS="$COMPFLAGS /openmp"']; - else - command = [command ' CFLAGS="\$CFLAGS -fopenmp"' ... - ' CXXFLAGS="\$CXXFLAGS -fopenmp"' ... - ' LDFLAGS="\$LDFLAGS -fopenmp"']; - end - case 'warnings' - if (ispc) - command = [command ' -v COMPFLAGS="$COMPFLAGS -W3"']; - else - command = [command ' -v' ... - ' CFLAGS="\$CFLAGS -Wall -Wextra"' ... - ' CXXFLAGS="\$CXXFLAGS -Wall -Wextra"']; - end - otherwise - error('Damkjer:makeVpsTree:NotImplementedFlag', ... - ['Flag recognized, but not supported (oops): ' flag]); - end % switch property - end - - varargin(1) = []; + flag = varargin{1}; + % If the property has been supplied in a shortened form, lengthen it + iFlag = find(strncmpi(flag, flags, length(flag))); + + if isempty(iFlag) + options = strcat(options, ' ', flag); + + % Pass through any arguments we don't recognize, alternatively escalate + % to an error. +% error('Damkjer:makeVpsTree:InvalidFlag', ['Invalid Flag: ' flag]); + elseif length(iFlag) > 1 + error('Damkjer:makeVpsTree:AmbiguousFlag', ... + ['Supplied shortened flag is ambiguous: ' flag]); + else + flag = flags{iFlag}; + + switch flag + case 'verbose' + command = strcat(command, ' -v'); + case 'debug' + command = strcat(command, ' -g'); + case 'parallel' + if (ispc) + command = strcat(command, ' COMPFLAGS="$COMPFLAGS /openmp"'); + else + command = strcat(command, ' CFLAGS="\$CFLAGS -fopenmp"',... + ' CXXFLAGS="\$CXXFLAGS -fopenmp"',... + ' LDFLAGS="\$LDFLAGS -fopenmp"'); + end + case 'warnings' + if (ispc) + command = strcat(command, ' -v COMPFLAGS="$COMPFLAGS -W3"'); + else + command = strcat(command, ' -v',... + ' CFLAGS="\$CFLAGS -Wall -Wextra"',... + ' CXXFLAGS="\$CXXFLAGS -Wall -Wextra"'); + end + otherwise + error('Damkjer:makeVpsTree:NotImplementedFlag', ... + ['Flag recognized, but not supported (oops): ' flag]); + end % switch property + end + + varargin(1) = []; end % while [path,~,~]=fileparts(mfilename('fullpath')); command = [command options ' -I' path '/../../.. -outdir ' path ' ' path '/']; + eval([command 'newVpTree.cpp']); eval([command 'deleteVpTree.cpp']); Index: Damkjer/Util/Widgets/timebar.m =================================================================== --- Damkjer/Util/Widgets/timebar.m (revision 0) +++ Damkjer/Util/Widgets/timebar.m (revision 6) @@ -1,58 +1,89 @@ -%%TIMEBAR Create a time-based waitbar. -% [HANDLE] = TIMEBAR(...) . +% timebar Create a time-based waitbar. % -% See also WAITBAR. +% File: +% timebar.m +% +% Description: +% This MATLAB function creates a waitbar dialog with an estimated time to +% completion based on linear extrapolation from number of processing steps +% completed, number of processing steps remaining and current elapsed time. +% +% Limitations: +% None. +% +% Synopsis: +% [handle] = timebar(step, steps, message, start, handle) +% +% Inputs: +% step - Current processing step. +% steps - Total steps in process. +% message - Base message to be included in waitbar text. +% start - Process start time. +% handle - Figure handle for existing waitbar. +% +% Outputs: +% handle - Figure handle for rendered timebar graphic. +% +% Toolbox requirements: +% None. +% +% Script requirements: +% None. +% +% Data requirements: +% None. +% +% References: +% None. +% +% See Also: +% waitbar +% -% Copyright 2012 Kristian Linn Damkjer +% Copyright (C) 2012 Kristian L. Damkjer. +% +% Software History: +% 2012-AUG-29 K. Damkjer +% Initial Coding. + function [handle] = timebar(step, steps, message, start, handle) -try - -step=double(step); -steps=double(steps); -elapsed=toc(start); -total=floor(elapsed/step*steps); -elapsed=floor(elapsed); -remain=total-elapsed; +step = double(step); +steps = double(steps); +elapsed = toc(start); +total = floor(elapsed/step*steps); +elapsed = floor(elapsed); +remain = total-elapsed; % Construct total string -if (total > 7200) - tmsg=[num2str(floor(total/3600)) 'h']; -elseif (total > 120) - tmsg=[num2str(floor(total/60)) 'm']; -else - tmsg=[num2str(total) 's']; -end +tmsg=timeString(total); % Construct elapsed string -if (elapsed > 7200) - emsg=[num2str(floor(elapsed/3600)) 'h']; -elseif (elapsed > 120) - emsg=[num2str(floor(elapsed/60)) 'm']; -else - emsg=[num2str(elapsed) 's']; -end +emsg=timeString(elapsed); % Construct total string -if (remain > 7200) - rmsg=[num2str(floor(remain/3600)) 'h']; -elseif (remain > 120) - rmsg=[num2str(floor(remain/60)) 'm']; +rmsg=timeString(remain); + +% Construct waitbar message +wbmsg = [message emsg '/' tmsg ' (' rmsg ' remaining)']; + +% Create a new waitbar if the handle is invalid or wasn't passed in. Otherwise +% render to the valid figure handle. +if (nargin > 4 && all(ishghandle(handle, 'figure'))) + waitbar(step/steps, handle, wbmsg); else - rmsg=[num2str(remain) 's']; -end - -if (nargin > 4 && all(ishghandle(handle, 'figure'))) - waitbar(step/steps,handle,... - [message emsg '/' tmsg ' (' rmsg ' remaining)']); -else - handle = waitbar(step/steps,... - [message emsg '/' tmsg ' (' rmsg ' remaining)']); -end - -catch me - disp(me.identifier); - disp(me.message); + handle = waitbar(step/steps, wbmsg); end end + + +function [message] = timeString(seconds) + if (seconds > 7200) + message=[num2str(floor(seconds/3600)) 'h']; + elseif (seconds > 120) + message=[num2str(floor(seconds/60)) 'm']; + else + message=[num2str(seconds) 's']; + end +end Index: Damkjer/applications/demo_dimcolor.m =================================================================== --- Damkjer/applications/demo_dimcolor.m (revision 0) +++ Damkjer/applications/demo_dimcolor.m (revision 6) @@ -6,4 +6,5 @@ function demo_dimcolor dataFile='Serpent Mound Model LAS Data.las'; +neighbors=20; disp('Building MEX files...'); @@ -11,30 +12,35 @@ makeMath('p'); makeVpTree('p'); -disp([' ...done: ' num2str(toc(buildStart)) ' s']); +disp([' ...done: ' timeFormat(toc(buildStart))]); -[path,~,~]=fileparts(mfilename('fullpath')); +[pathstr,~,~] = fileparts(mfilename('fullpath')); +[pathstr,~,~] = fileparts(pathstr); -srcPath = [path '/../data/source']; -outPath = [path '/../data/results']; +srcPath = fullfile(pathstr, 'data', 'source'); +outPath = fullfile(pathstr, 'data', 'results'); disp('Loading Data...'); loadStart=tic; -sourceData=LasFile([srcPath '/' dataFile]); -disp([' ...done: ' num2str(toc(loadStart)) ' s']); +sourceData=LasFile(fullfile(srcPath, dataFile)); +disp([' ...done: read ' num2str(size(sourceData.points,2)) ... + ' points in ' timeFormat(toc(loadStart))]); disp('Building Database...'); dbStart=tic; database=VpTree(sourceData.points); -disp([' ...done: ' num2str(toc(dbStart)) ' s']); +disp([' ...done: ' timeFormat(toc(dbStart))]); disp('Colorizing by Inherent Dimensionality...'); colorStart=tic; -sourceData.color=dimcolor(sourceData.points,database,'n',27); -disp([' ...done: ' num2str(toc(colorStart)) ' s']); +sourceData.color=dimcolor(sourceData.points,database,'n',neighbors); +disp([' ...done: evaluated ' num2str(size(sourceData.points,2)) ... + ' points in ' num2str(neighbors) '-point neighborhoods in ' ... + timeFormat(toc(colorStart))]); -outFile=[outPath '/dimcolor.las']; +outFile=fullfile(outPath, 'dimcolor.las'); disp(['Saving Results to: ' outFile '...']); writeStart=tic; sourceData.saveTo(outFile); -disp([' ...done: ' num2str(toc(writeStart)) ' s']); +disp([' ...done: wrote ' num2str(size(sourceData.points,2)) ... + ' points in ' timeFormat(toc(writeStart))]); end Index: Damkjer/applications/demo_spatial_analyze.m =================================================================== --- Damkjer/applications/demo_spatial_analyze.m (revision 6) +++ Damkjer/applications/demo_spatial_analyze.m (revision 6) @@ -0,0 +1,50 @@ +%% +% DEMO_SPATIAL_ANALYZE +% +% Demonstrates the spatial analyzer function capability. +% +function demo_spatial_analyze +%dataFile='Serpent Mound Model LAS Data.las'; +dataFile='1815-3_Armstrong-Enderby_Sample1.las'; +neighbors=20; +%counts=10:20; + +disp('Building MEX files...'); +buildStart=tic; +makeMath('p'); +makeVpTree('p'); +disp([' ...done: ' num2str(toc(buildStart)) 's']); + +[path,~,~]=fileparts(mfilename('fullpath')); + +srcPath = [path '/../data/source']; +outPath = [path '/../data/results']; + +disp('Loading Data...'); +loadStart=tic; +sourceData=LasFile([srcPath '/' dataFile]); +disp([' ...done: read ' num2str(size(sourceData.points,2)) ' points in ' ... + num2str(toc(loadStart)) 's']); + +disp('Building Database...'); +dbStart=tic; +database=VpTree(sourceData.points); +disp([' ...done: ' num2str(toc(dbStart)) 's']); + +disp('Performing Spatial Analysis...'); +resultsStart=tic; +results=spanalyze(sourceData.points,database,'n',neighbors); +%results=spanalyze(sourceData.points,database,'c',counts); +disp([' ...done: evaluated ' num2str(size(sourceData.points,2)) ' points in ' ... + num2str(neighbors) '-point neighborhoods in ' ... + num2str(toc(resultsStart)) 's']); + +sourceData.color=[uint16(double(intmax('uint16'))*results.dimensionality)]; + +outFile=[outPath '/dimcolor.las']; +disp(['Saving Results to: ' outFile '...']); +writeStart=tic; +sourceData.saveTo(outFile); +disp([' ...done: wrote ' num2str(size(sourceData.points,2)) ' points in ' ... + num2str(toc(writeStart)) 's']); +end