1 {
2 lasvectorialreader.pas
3 
4 Reads geospatial data encoded in the ASPRS LASer (LAS) file format version 1.3
5 
6 LAS file format specification obtained from:
7 
8 http://www.asprs.org/a/society/committees/standards/lidar_exchange_format.html
9 LAS 1.3 r11
10 
11 All data is in little-endian format.
12 
13 AUTHORS: Felipe Monteiro de Carvalho
14 
15 License: The same modified LGPL as the Free Pascal RTL
16          See the file COPYING.modifiedLGPL for more details
17 }
18 unit lasvectorialreader;
19 
20 {$ifdef fpc}
21   {$mode delphi}
22 {$endif}
23 
24 {$define FPVECTORIALDEBUG_LAS}
25 
26 interface
27 
28 uses
29   Classes, SysUtils, dateutils,
30   fpcanvas, fpimage,
31   //avisozlib,
32   fpvectorial;
33 
34 type
35   // LAS data types introduced in LAS 1.0
36   laschar = AnsiChar; // or ShortInt
37   lasuchar = Byte;
38   lasshort = Smallint;
39   lasushort = Word;
40   laslong = Integer;
41   lasulong = Cardinal;
42   lasdouble = double;
43   //
44   laslonglong = Int64;
45   lasulonglong = QWord;
46 
47   // PUBLIC HEADER BLOCK version 1.0
48   TLASPublicHeaderBlock_1_0 = packed record
49     FileSignatureLASF: array[0..3] of laschar;
50     FileSourceID: lasushort;                    // Reserved in LAS 1.0
51     GlobalEncoding: lasushort;                  // Reserved in LAS 1.0
52     ProjectIDGUIDdata1: lasulong;               // Optional
53     ProjectIDGUIDdata2: lasushort;              // Optional
54     ProjectIDGUIDdata3: lasushort;              // Optional
55     ProjectIDGUIDdata4: array[0..7] of lasuchar;// Optional
56     VersionMajor: lasuchar;
57     VersionMinor: lasuchar;
58     SystemIdentifier: array[0..31] of laschar;
59     GeneratingSoftware: array[0..31] of laschar;
60     FileCreationDayofYear: lasushort;           // Name in LAS 1.0 -> FlightDateJulian, but same meaning
61     FileCreationYear: lasushort;                // Name in LAS 1.0 -> Year, but same meaning
62     HeaderSize: lasushort;                      // Name in LAS 1.0 -> OffsetToData
63     OffsetToPointData: lasulong;
64     NumberofVariableLengthRecords: lasulong;
65     PointDataFormatID: lasuchar; // (0-99 for spec)
66     PointDataRecordLength: lasushort;
67     Numberofpointrecords: lasulong;
68     Numberofpointsbyreturn: array[0..4] of lasulong;
69     Xscalefactor: lasdouble;
70     Yscalefactor: lasdouble;
71     Zscalefactor: lasdouble;
72 
73     Xoffset: lasdouble;
74     Yoffset: lasdouble;
75     Zoffset: lasdouble;
76     MaxX: lasdouble;
77     MinX: lasdouble;
78     MaxY: lasdouble;
79     MinY: lasdouble;
80     MaxZ: lasdouble;
81     MinZ: lasdouble;
82   end;
83 
84   // PUBLIC HEADER BLOCK Extension in Version 1.3
85   TLASPublicHeaderBlock_1_3_Extension = packed record
86     StartofWaveformDataPacket: lasulonglong;
87   end;
88 
89   TLASVariableLengthRecord = packed record
90     RecordSignatureAABB: lasushort;
91     UserID: array[0..15] of laschar;
92     RecordID: lasushort;
93     RecordLengthAfterHeader: lasushort;
94     Description: array[0..31] of laschar;
95   end;
96 
97   // Points are split in the following atoms for compression:
98   // POINT10, GPSTIME10, RGB12, WAVEPACKET13, and BYTE
99   // for the LAZ format
100 
101   // Contains POINT10
102   TLASPointDataRecordFormat0 = packed record
103     X: laslong;
104     Y: laslong;
105     Z: laslong;
106     Intensity: lasushort;
107     Flags: Byte;
108     Classification: lasuchar;
109     ScanAngleRank: lasuchar; // (-90 to +90) - is the left side
110     FileMarker: lasuchar;
111     UserBitField: lasushort;
112   end;
113 
114   // Contains POINT10 + GPSTIME10
115   TLASPointDataRecordFormat1 = packed record
116     X: laslong;
117     Y: laslong;
118     Z: laslong;
119     Intensity: lasushort;
120     Flags: Byte;
121     Classification: lasuchar;
122     ScanAngleRank: lasuchar; // (-90 to +90) - is the left side
123     FileMarker: lasuchar;
124     UserBitField: lasushort;
125     GPSTime: lasdouble;
126   end;
127 
128   { TvLASVectorialReader }
129 
130   TvLASVectorialReader = class(TvCustomVectorialReader)
131   protected
132     // Stream position information
133     InitialPos, PositionAfterPublicHeader: Int64;
134     {$ifdef FPVECTORIALDEBUG_LAS}
135     procedure DebugOutPublicHeaderBlock();
136     {$endif}
137     procedure ReadVariableLengthRecords(AStream: TStream);
138     procedure DoProgress(AProgress: Byte; AData: TvVectorialDocument);
ReadLAZPoint0null139     function ReadLAZPoint0(AStream: TStream): TLASPointDataRecordFormat0;
140   public
141     // Public Header
142     PublicHeaderBlock_1_0: TLASPublicHeaderBlock_1_0;
143     PublicHeaderBlock_1_3_Extension: TLASPublicHeaderBlock_1_3_Extension;
144     // Variable Length Records
145     VariableLengthRecords: array of TLASVariableLengthRecord;
146     // Point Data
147     PointsFormat0: array of TLASPointDataRecordFormat0;
148     PointsFormat1: array of TLASPointDataRecordFormat1;
149     { General reading methods }
150     procedure ReadFromStream(AStream: TStream; AData: TvVectorialDocument); override;
151   end;
152 
153   { TvLASVectorialWriter }
154 
155   TvLASVectorialWriter = class(TvCustomVectorialWriter)
156   private
157     // Stream position information
158     InitialPos, PositionAfterPublicHeader: Int64;
159   public
160     // Public Header
161     PublicHeaderBlock_1_0: TLASPublicHeaderBlock_1_0;
162     PublicHeaderBlock_1_3_Extension: TLASPublicHeaderBlock_1_3_Extension;
163     // Variable Length Records
164     VariableLengthRecords: array of TLASVariableLengthRecord;
165     // Point Data
166     PointsFormat0: array of TLASPointDataRecordFormat0;
167     PointsFormat1: array of TLASPointDataRecordFormat1;
168     { General reading methods }
169     procedure WriteToStream(AStream: TStream; AData: TvVectorialDocument); override;
170   end;
171 
172 implementation
173 
174 { TvLASVectorialWriter }
175 
176 procedure TvLASVectorialWriter.WriteToStream(AStream: TStream;
177   AData: TvVectorialDocument);
178 var
179   lPage: TvVectorialPage;
180   lRecord0: TLASPointDataRecordFormat0;
181   lRecord1: TLASPointDataRecordFormat1;
182   lPoint: TvPoint;
183   lColor: TFPColor;
184   lCreationDate: TDateTime;
185   lEntity: TvEntity;
186   i: Integer;
187 begin
188   // Get the first page
189   lPage := AData.GetPageAsVectorial(0);
190   lCreationDate := Now;
191 
192   // Write our LAS 1.0 header
193   FillChar(PublicHeaderBlock_1_0, SizeOf(PublicHeaderBlock_1_0), #0);
194   PublicHeaderBlock_1_0.FileSignatureLASF := 'LASF';
195   PublicHeaderBlock_1_0.FileSourceID := 0;
196   PublicHeaderBlock_1_0.GlobalEncoding := 0;
197   PublicHeaderBlock_1_0.ProjectIDGUIDdata1 := 0;
198   PublicHeaderBlock_1_0.ProjectIDGUIDdata2 := 0;
199   PublicHeaderBlock_1_0.ProjectIDGUIDdata3 := 0;
200   // PublicHeaderBlock_1_0.ProjectIDGUIDdata4 all zero
201   PublicHeaderBlock_1_0.VersionMajor := 1;
202   PublicHeaderBlock_1_0.VersionMinor := 0;
203   PublicHeaderBlock_1_0.SystemIdentifier := '';
204   PublicHeaderBlock_1_0.GeneratingSoftware := 'FPSpreadsheet';
205   PublicHeaderBlock_1_0.FileCreationDayofYear := DayOfTheYear(lCreationDate);
206   PublicHeaderBlock_1_0.FileCreationYear := YearOf(lCreationDate);
207   PublicHeaderBlock_1_0.HeaderSize := SizeOf(PublicHeaderBlock_1_0);
208   PublicHeaderBlock_1_0.OffsetToPointData := SizeOf(PublicHeaderBlock_1_0);
209   PublicHeaderBlock_1_0.NumberofVariableLengthRecords := 0;
210   PublicHeaderBlock_1_0.PointDataFormatID := 1;
211   PublicHeaderBlock_1_0.PointDataRecordLength := $1C;
212   PublicHeaderBlock_1_0.Numberofpointrecords := lPage.GetEntitiesCount;
213   PublicHeaderBlock_1_0.Numberofpointsbyreturn[0] := 0;
214   PublicHeaderBlock_1_0.Numberofpointsbyreturn[1] := 0;
215   PublicHeaderBlock_1_0.Numberofpointsbyreturn[2] := 0;
216   PublicHeaderBlock_1_0.Numberofpointsbyreturn[3] := 0;
217   PublicHeaderBlock_1_0.Numberofpointsbyreturn[4] := 0;
218   PublicHeaderBlock_1_0.Xscalefactor := 1;
219   PublicHeaderBlock_1_0.Yscalefactor := 1;
220   PublicHeaderBlock_1_0.Zscalefactor := 1;
221 
222   PublicHeaderBlock_1_0.Xoffset := 0;
223   PublicHeaderBlock_1_0.Yoffset := 0;
224   PublicHeaderBlock_1_0.Zoffset := 0;
225   PublicHeaderBlock_1_0.MaxX := lPage.MaxX;
226   PublicHeaderBlock_1_0.MinX := lPage.MinX;
227   PublicHeaderBlock_1_0.MaxY := lPage.MaxY;
228   PublicHeaderBlock_1_0.MinY := lPage.MinY;
229   PublicHeaderBlock_1_0.MaxZ := lPage.MaxZ;
230   PublicHeaderBlock_1_0.MinZ := lPage.MinZ;
231   AStream.Write(PublicHeaderBlock_1_0, SizeOf(TLASPublicHeaderBlock_1_0));
232 
233   // Write the variable length records
234   // none currently
235 
236   // Write the point data
237   for i := 0 to lPage.GetEntitiesCount()-1 do
238   begin
239     lEntity := lPage.GetEntity(i);
240     if not (lEntity is TvPoint) then Continue;
241     lPoint := lEntity as TvPoint;
242 
243     FillChar(lRecord1, SizeOf(TLASPointDataRecordFormat1), #0);
244     lRecord1.X := Round(lEntity.X);
245     lRecord1.Y := Round(lEntity.Y);
246     lRecord1.Z := Round(lEntity.Z);
247 
248     // Convert the colors into LIDAR Point Classes
249     lColor := lPoint.Pen.Color;
250     // 2 Ground
251     if lColor = colMaroon then
252       lRecord1.Classification := 2
253     // 3 Low Vegetation
254     else if lColor = colGreen then
255       lRecord1.Classification := 3
256     // 4 Medium Vegetation
257     //4: lColor := colGreen;
258     // 5 High Vegetation
259     else if lColor = colDkGreen then
260       lRecord1.Classification := 5
261     // 6 Building
262     else if lColor = colGray then
263       lRecord1.Classification := 6
264     // 7 Low Point (noise)
265     // 8 Model Key-point (mass point)
266     // 9 Water
267     else if lColor = colBlue then
268       lRecord1.Classification := 9;
269 
270     AStream.Write(lRecord1, SizeOf(TLASPointDataRecordFormat1));
271   end;
272 end;
273 
274 {$ifdef FPVECTORIALDEBUG_LAS}
275 procedure TvLASVectorialReader.DebugOutPublicHeaderBlock;
276 begin
277   WriteLn(Format('FileSignatureLASF = %s = %x %x %x %x',
278     [string(PublicHeaderBlock_1_0.FileSignatureLASF),
279      Ord(PublicHeaderBlock_1_0.FileSignatureLASF[0]),
280      Ord(PublicHeaderBlock_1_0.FileSignatureLASF[1]),
281      Ord(PublicHeaderBlock_1_0.FileSignatureLASF[2]),
282      Ord(PublicHeaderBlock_1_0.FileSignatureLASF[3])]));
283   WriteLn(Format('FileSourceID = $%x', [PublicHeaderBlock_1_0.FileSourceID]));
284   WriteLn(Format('GlobalEncoding = $%x', [PublicHeaderBlock_1_0.GlobalEncoding]));
285   WriteLn(Format('ProjectIDGUIDdata1 = $%x', [PublicHeaderBlock_1_0.ProjectIDGUIDdata1]));
286   WriteLn(Format('ProjectIDGUIDdata2 = $%x', [PublicHeaderBlock_1_0.ProjectIDGUIDdata2]));
287   WriteLn(Format('ProjectIDGUIDdata3 = $%x', [PublicHeaderBlock_1_0.ProjectIDGUIDdata3]));
288 //  WriteLn(Format('ProjectIDGUIDdata4 = %x', [ProjectIDGUIDdata2]));
289   WriteLn(Format('VersionMajor = %d', [PublicHeaderBlock_1_0.VersionMajor]));
290   WriteLn(Format('VersionMinor = %d', [PublicHeaderBlock_1_0.VersionMinor]));
291   WriteLn(Format('SystemIdentifier = %s', [PublicHeaderBlock_1_0.SystemIdentifier]));
292   WriteLn(Format('GeneratingSoftware = %s', [PublicHeaderBlock_1_0.GeneratingSoftware]));
293   WriteLn(Format('FileCreationDayofYear = %d', [PublicHeaderBlock_1_0.FileCreationDayofYear]));
294   WriteLn(Format('FileCreationYear = %d', [PublicHeaderBlock_1_0.FileCreationYear]));
295   WriteLn(Format('HeaderSize = $%x', [PublicHeaderBlock_1_0.HeaderSize]));
296   WriteLn(Format('OffsetToPointData = $%x', [PublicHeaderBlock_1_0.OffsetToPointData]));
297   WriteLn(Format('NumberofVariableLengthRecords = $%x', [PublicHeaderBlock_1_0.NumberofVariableLengthRecords]));
298   WriteLn(Format('PointDataFormatID = $%x', [PublicHeaderBlock_1_0.PointDataFormatID]));
299   WriteLn(Format('PointDataRecordLength = $%x', [PublicHeaderBlock_1_0.PointDataRecordLength]));
300   WriteLn(Format('Numberofpointrecords = $%x', [PublicHeaderBlock_1_0.Numberofpointrecords]));
301   WriteLn(Format('Numberofpointsbyreturn = %x %x %x %x %x',
302     [PublicHeaderBlock_1_0.Numberofpointsbyreturn[0],
303      PublicHeaderBlock_1_0.Numberofpointsbyreturn[1],
304      PublicHeaderBlock_1_0.Numberofpointsbyreturn[2],
305      PublicHeaderBlock_1_0.Numberofpointsbyreturn[3],
306      PublicHeaderBlock_1_0.Numberofpointsbyreturn[4]
307      ]));
308   WriteLn(Format('Xscalefactor = %f', [PublicHeaderBlock_1_0.Xscalefactor]));
309   WriteLn(Format('Yscalefactor = %f', [PublicHeaderBlock_1_0.Yscalefactor]));
310   WriteLn(Format('Zscalefactor = %f', [PublicHeaderBlock_1_0.Zscalefactor]));
311   WriteLn(Format('Xoffset = %f', [PublicHeaderBlock_1_0.Xoffset]));
312   WriteLn(Format('Yoffset = %f', [PublicHeaderBlock_1_0.Yoffset]));
313   WriteLn(Format('Zoffset = %f', [PublicHeaderBlock_1_0.Zoffset]));
314   WriteLn(Format('MaxX = %f', [PublicHeaderBlock_1_0.MaxX]));
315   WriteLn(Format('MinX = %f', [PublicHeaderBlock_1_0.MinX]));
316   WriteLn(Format('MaxY = %f', [PublicHeaderBlock_1_0.MaxY]));
317   WriteLn(Format('MinY = %f', [PublicHeaderBlock_1_0.MinY]));
318   WriteLn(Format('MaxZ = %f', [PublicHeaderBlock_1_0.MaxZ]));
319   WriteLn(Format('MinZ = %f', [PublicHeaderBlock_1_0.MinZ]));
320   WriteLn('');
321   WriteLn(Format('LAS 1.0 header size = $%x', [SizeOf(TLASPublicHeaderBlock_1_0)]));
322   WriteLn('');
323 end;
324 {$endif}
325 
326 procedure TvLASVectorialReader.ReadVariableLengthRecords(AStream: TStream);
327 var
328   i: Integer;
329   NextPosition: Int64;
330 begin
331   NextPosition := PositionAfterPublicHeader;
332 
333   SetLength(VariableLengthRecords, PublicHeaderBlock_1_0.NumberofVariableLengthRecords);
334   for i := 0 to PublicHeaderBlock_1_0.NumberofVariableLengthRecords - 1 do
335   begin
336     AStream.Position := NextPosition;
337     {$ifdef FPVECTORIALDEBUG_LAS}
338     WriteLn(Format('Variable Length Record #%d Position = $%x', [i, NextPosition]));
339     WriteLn('');
340     {$endif}
341     AStream.Read(VariableLengthRecords[i], SizeOf(TLASVariableLengthRecord));
342     NextPosition := AStream.Position+VariableLengthRecords[i].RecordLengthAfterHeader;
343     {$ifdef FPVECTORIALDEBUG_LAS}
344     WriteLn(Format('RecordSignatureAABB = $%x', [VariableLengthRecords[i].RecordSignatureAABB]));
345     WriteLn(Format('UserID = %s', [VariableLengthRecords[i].UserID]));
346     WriteLn(Format('RecordID = $%x', [VariableLengthRecords[i].RecordID]));
347     WriteLn(Format('RecordLengthAfterHeader = $%x', [VariableLengthRecords[i].RecordLengthAfterHeader]));
348     WriteLn(Format('Description = %s', [VariableLengthRecords[i].Description]));
349     WriteLn('');
350     {$endif}
351   end;
352 end;
353 
354 procedure TvLASVectorialReader.DoProgress(AProgress: Byte; AData: TvVectorialDocument);
355 begin
356   if @AData.OnProgress <> nil then AData.OnProgress(AProgress);
357 end;
358 
TvLASVectorialReader.ReadLAZPoint0null359 function TvLASVectorialReader.ReadLAZPoint0(AStream: TStream
360   ): TLASPointDataRecordFormat0;
361 begin
362 
363 end;
364 
365 procedure TvLASVectorialReader.ReadFromStream(AStream: TStream;
366   AData: TvVectorialDocument);
367 var
368   lPage: TvVectorialPage;
369   lRecord0: TLASPointDataRecordFormat0;
370   lRecord1: TLASPointDataRecordFormat1;
371   lFirstPoint: Boolean = True;
372   lPoint: TvPoint;
373   lClassification: Integer = -1;
374   lColor: TFPColor;
375   lPointsCounter: Integer = 0;
376 begin
377   // Clear and add the first page
378   AData.Clear;
379   lPage := AData.AddPage();
380 
381   // First read the header like if it was for LAS 1.0,
382   // this will tell us the real version so then we read it again
383   InitialPos := AStream.Position;
384   AStream.Read(PublicHeaderBlock_1_0, SizeOf(TLASPublicHeaderBlock_1_0));
385   {$ifdef FPVECTORIALDEBUG_LAS}
386   DebugOutPublicHeaderBlock();
387   {$endif}
388 
389   // First check the signature
390   if PublicHeaderBlock_1_0.FileSignatureLASF <> 'LASF' then
391     raise Exception.Create('[TvLASVectorialReader.ReadFromStream] Invalid file signoture while reading LAS file');
392 
393   lPage.MinX := PublicHeaderBlock_1_0.MinX;
394   lPage.MinY := PublicHeaderBlock_1_0.MinY;
395   lPage.MinZ := PublicHeaderBlock_1_0.MinZ;
396   lPage.MaxX := PublicHeaderBlock_1_0.MaxX;
397   lPage.MaxY := PublicHeaderBlock_1_0.MaxY;
398   lPage.MaxZ := PublicHeaderBlock_1_0.MaxZ;
399 
400   // In LAS 1.3+ read the header extension
401   // ToDo
402 
403   PositionAfterPublicHeader := AStream.Position;
404 
405   // Read the variable length records
406   ReadVariableLengthRecords(AStream);
407 
408   // Read the point data
409   AStream.Position := InitialPos + PublicHeaderBlock_1_0.OffsetToPointData;
410   while AStream.Position < AStream.Size do
411   begin
412     // Send a progress event every 1k points
413     Inc(lPointsCounter);
414     if lPointsCounter mod 1000 = 0 then
415       DoProgress(Round(AStream.Position * 100 / AStream.Size), AData);
416 
417     // hack to cut las files: if lPage.GetEntitiesCount = 100000 then Exit;
418     case PublicHeaderBlock_1_0.PointDataFormatID of
419     0:
420     begin
421       AStream.ReadBuffer(lRecord0, SizeOf(TLASPointDataRecordFormat0));
422       lClassification := lRecord0.Classification;
423       lPoint := lPage.AddPoint(lRecord0.X, lRecord0.Y, lRecord0.Z);
424     end;
425     1:
426     begin
427       AStream.ReadBuffer(lRecord1, SizeOf(TLASPointDataRecordFormat1));
428       lClassification := lRecord1.Classification;
429       lPoint := lPage.AddPoint(lRecord1.X, lRecord1.Y, lRecord1.Z);
430     end;
431     130:
432     begin
433       lRecord0 := ReadLAZPoint0(AStream);
434       lClassification := lRecord0.Classification;
435       lPoint := lPage.AddPoint(lRecord0.X, lRecord0.Y, lRecord0.Z);
436     end;
437     else
438       raise Exception.Create('[TvLASVectorialReader.ReadFromStream] Error reading LAS point: Unknown point type number='
439        + IntToStr(PublicHeaderBlock_1_0.PointDataFormatID));
440     end;
441 
442     // Correct the min and max
443     if lFirstPoint then
444     begin
445       lPage.MinX := lPoint.X;
446       lPage.MinY := lPoint.Y;
447       lPage.MinZ := lPoint.Z;
448       lPage.MaxX := lPoint.X;
449       lPage.MaxY := lPoint.Y;
450       lPage.MaxZ := lPoint.Z;
451 
452       lFirstPoint := False;
453     end
454     else
455     begin
456       if lPage.MinX > lPoint.X then lPage.MinX := lPoint.X;
457       if lPage.MinY > lPoint.Y then lPage.MinY := lPoint.Y;
458       if lPage.MinZ > lPoint.Z then lPage.MinZ := lPoint.Z;
459       if lPage.MaxX < lPoint.X then lPage.MaxX := lPoint.X;
460       if lPage.MaxY < lPoint.Y then lPage.MaxY := lPoint.Y;
461       if lPage.MaxZ < lPoint.Z then lPage.MaxZ := lPoint.Z;
462     end;
463 
464     // Set a color if desired
465     case lClassification of
466     // Table 4.9 - ASPRS Standard LIDAR Point Classes
467     // Classification Value
468     // 0 Created, never classified
469     // 1 Unclassified1
470     // 2 Ground
471     2: lColor := colMaroon;
472     // 3 Low Vegetation
473     3: lColor := colGreen;
474     // 4 Medium Vegetation
475     4: lColor := colGreen;
476     // 5 High Vegetation
477     5: lColor := colDkGreen;
478     // 6 Building
479     6: lColor := colGray;
480     // 7 Low Point (noise)
481     // 8 Model Key-point (mass point)
482     // 9 Water
483     9: lColor := colBlue;
484     else
485       lColor := colBlack;
486     end;
487 
488     lPoint.Pen.Color := lColor;
489   end;
490 end;
491 
492 initialization
493 
494   RegisterVectorialReader(TvLASVectorialReader, vfLAS);
495   RegisterVectorialWriter(TvLASVectorialWriter, vfLAS);
496   RegisterVectorialReader(TvLASVectorialReader, vfLAZ);
497 
498 end.
499 
500