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