1 /******************************************************************************
2 * Copyright (c) 2011, Michael P. Gerlek (mpg@flaxen.com)
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following
8 * conditions are met:
9 *
10 * * Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * * Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided
15 * with the distribution.
16 * * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
17 * names of its contributors may be used to endorse or promote
18 * products derived from this software without specific prior
19 * written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
24 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
25 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
27 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
28 * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
29 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
30 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
31 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
32 * OF SUCH DAMAGE.
33 ****************************************************************************/
34
35 #include <pdal/compression/LazPerfVlrCompression.hpp>
36
37 #include "LasHeader.hpp"
38 #include "LasReader.hpp"
39 #include "LasUtils.hpp"
40
41 #include <sstream>
42 #include <string.h>
43
44 #include <pdal/pdal_features.hpp>
45 #include <pdal/Metadata.hpp>
46 #include <pdal/PointView.hpp>
47 #include <pdal/QuickInfo.hpp>
48 #include <pdal/util/Extractor.hpp>
49 #include <pdal/util/FileUtils.hpp>
50 #include <pdal/util/IStream.hpp>
51 #include <pdal/util/ProgramArgs.hpp>
52
53 #include "LasHeader.hpp"
54 #include "LasVLR.hpp"
55
56 #ifdef PDAL_HAVE_LASZIP
57 #include <laszip/laszip_api.h>
58 #else
59 using laszip_POINTER = void *;
60 using laszip_point_struct = void *;
61 struct laszip_point;
62 #endif
63
64 namespace pdal
65 {
66
67 namespace
68 {
69
70 struct invalid_stream : public std::runtime_error
71 {
invalid_streampdal::__anon618467ac0111::invalid_stream72 invalid_stream(const std::string& msg) : std::runtime_error(msg)
73 {}
74 };
75
76 } // unnamed namespace
77
78 struct LasReader::Args
79 {
80 public:
81 StringList extraDimSpec;
82 std::string compression;
83 bool useEbVlr;
84 StringList ignoreVLROption;
85 bool fixNames;
86 PointId start;
87 bool nosrs;
88 };
89
90 struct LasReader::Private
91 {
92 typedef std::vector<LasUtils::IgnoreVLR> IgnoreVLRList;
93
94 LasHeader header;
95 laszip_POINTER laszip;
96 laszip_point_struct *laszipPoint;
97 LazPerfVlrDecompressor *decompressor;
98 std::vector<char> decompressorBuf;
99 point_count_t index;
100 IgnoreVLRList ignoreVLRs;
101 std::vector<ExtraDim> extraDims;
102
Privatepdal::LasReader::Private103 Private() : decompressor(nullptr), index(0)
104 {}
105 };
106
LasReader()107 LasReader::LasReader() : m_args(new Args), m_p(new Private)
108 {}
109
110
~LasReader()111 LasReader::~LasReader()
112 {
113 #ifdef PDAL_HAVE_LAZPERF
114 delete m_p->decompressor;
115 #endif
116 }
117
118
addArgs(ProgramArgs & args)119 void LasReader::addArgs(ProgramArgs& args)
120 {
121 args.add("extra_dims", "Dimensions to assign to extra byte data",
122 m_args->extraDimSpec);
123 args.add("compression", "Decompressor to use", m_args->compression, "EITHER");
124 args.add("use_eb_vlr", "Use extra bytes VLR for 1.0 - 1.3 files", m_args->useEbVlr);
125 args.add("ignore_vlr", "VLR userid/recordid to ignore", m_args->ignoreVLROption);
126 args.add("start", "Point at which reading should start (0-indexed).", m_args->start);
127 args.add("fix_dims", "Make invalid dimension names valid by changing "
128 "invalid characters to '_'", m_args->fixNames, true);
129 args.add("nosrs", "Skip reading/processing file SRS", m_args->nosrs);
130 }
131
132
133 static StaticPluginInfo const s_info {
134 "readers.las",
135 "ASPRS LAS 1.0 - 1.4 read support. LASzip support is also \n" \
136 "enabled through this driver if LASzip was found during \n" \
137 "compilation.",
138 "http://pdal.io/stages/readers.las.html",
139 { "las", "laz" }
140 };
141
CREATE_STATIC_STAGE(LasReader,s_info)142 CREATE_STATIC_STAGE(LasReader, s_info)
143
144 std::string LasReader::getName() const { return s_info.name; }
145
header() const146 const LasHeader& LasReader::header() const
147 {
148 return m_p->header;
149 }
150
getNumPoints() const151 point_count_t LasReader::getNumPoints() const
152 {
153 return m_p->header.pointCount() - m_args->start;
154 }
155
initialize(PointTableRef table)156 void LasReader::initialize(PointTableRef table)
157 {
158 initializeLocal(table, m_metadata);
159 }
160
inspect()161 QuickInfo LasReader::inspect()
162 {
163 QuickInfo qi;
164 std::unique_ptr<PointLayout> layout(new PointLayout());
165
166 RowPointTable table;
167 initialize(table);
168 addDimensions(layout.get());
169
170 Dimension::IdList dims = layout->dims();
171 for (auto di = dims.begin(); di != dims.end(); ++di)
172 qi.m_dimNames.push_back(layout->dimName(*di));
173 if (!Utils::numericCast(m_p->header.pointCount(), qi.m_pointCount))
174 qi.m_pointCount = (std::numeric_limits<point_count_t>::max)();
175 qi.m_bounds = m_p->header.getBounds();
176 qi.m_srs = getSpatialReference();
177 qi.m_valid = true;
178
179 done(table);
180
181 return qi;
182 }
183
184
handleCompressionOption()185 void LasReader::handleCompressionOption()
186 {
187 std::string compression = Utils::toupper(m_args->compression);
188 #if defined(PDAL_HAVE_LAZPERF) && defined(PDAL_HAVE_LASZIP)
189 if (compression == "EITHER")
190 compression = "LASZIP";
191 #endif
192 #if !defined(PDAL_HAVE_LAZPERF) && defined(PDAL_HAVE_LASZIP)
193 if (compression == "EITHER")
194 compression = "LASZIP";
195 if (compression == "LAZPERF")
196 throwError("Can't decompress with LAZperf. PDAL not built "
197 "with LAZperf.");
198 #endif
199 #if defined(PDAL_HAVE_LAZPERF) && !defined(PDAL_HAVE_LASZIP)
200 if (compression == "EITHER")
201 compression = "LAZPERF";
202 if (compression == "LASZIP")
203 throwError("Can't decompress with LASzip. PDAL not built "
204 "with LASzip.");
205 #endif
206
207 #if defined(PDAL_HAVE_LAZPERF) || defined(PDAL_HAVE_LASZIP)
208 if (compression != "LAZPERF" && compression != "LASZIP")
209 throwError("Invalid value for option for compression: '" +
210 m_args->compression + "'. Value values are 'lazperf' and 'laszip'.");
211 #endif
212
213 // Set case-corrected value.
214 m_args->compression = compression;
215 }
216
createStream()217 void LasReader::createStream()
218 {
219 if (m_streamIf)
220 std::cerr << "Attempt to create stream twice!\n";
221 m_streamIf.reset(new LasStreamIf(m_filename));
222 if (!m_streamIf->m_istream)
223 {
224 std::ostringstream oss;
225 oss << "Unable to open stream for '"
226 << m_filename <<"' with error '" << strerror(errno) <<"'";
227 throw pdal_error(oss.str());
228 }
229 }
230
231
initializeLocal(PointTableRef table,MetadataNode & m)232 void LasReader::initializeLocal(PointTableRef table, MetadataNode& m)
233 {
234 try
235 {
236 m_p->extraDims = LasUtils::parse(m_args->extraDimSpec, false);
237 }
238 catch (const LasUtils::error& err)
239 {
240 throwError(err.what());
241 }
242
243 try
244 {
245 m_p->ignoreVLRs = LasUtils::parseIgnoreVLRs(m_args->ignoreVLROption);
246 }
247 catch (const LasUtils::error& err)
248 {
249 throwError(err.what());
250 }
251
252 m_p->header.initialize(log(), Utils::fileSize(m_filename), m_args->nosrs);
253 createStream();
254 std::istream *stream(m_streamIf->m_istream);
255
256 stream->seekg(0);
257 ILeStream in(stream);
258 try
259 {
260 // This also reads the extended VLRs at the end of the data.
261 in >> m_p->header;
262 }
263 catch (const LasHeader::error& e)
264 {
265 throwError(e.what());
266 }
267
268 for (auto i: m_p->ignoreVLRs)
269 {
270 if (i.m_recordId)
271 m_p->header.removeVLR(i.m_userId, i.m_recordId);
272 else
273 m_p->header.removeVLR(i.m_userId);
274 }
275
276 if (m_args->start > m_p->header.pointCount())
277 throwError("'start' value of " + std::to_string(m_args->start) + " is too large. "
278 "File contains " + std::to_string(m_p->header.pointCount()) + " points.");
279 if (m_p->header.compressed())
280 handleCompressionOption();
281 #ifdef PDAL_HAVE_LASZIP
282 m_p->laszip = nullptr;
283 #endif
284
285 if (!m_p->header.pointFormatSupported())
286 throwError("Unsupported LAS input point format: " +
287 Utils::toString((int)m_p->header.pointFormat()) + ".");
288
289 if (m_p->header.versionAtLeast(1, 4) || m_args->useEbVlr)
290 readExtraBytesVlr();
291 setSrs(m);
292 MetadataNode forward = table.privateMetadata("lasforward");
293 extractHeaderMetadata(forward, m);
294 extractVlrMetadata(forward, m);
295
296 m_streamIf.reset();
297 }
298
299
handleLaszip(int result)300 void LasReader::handleLaszip(int result)
301 {
302 #ifdef PDAL_HAVE_LASZIP
303 if (result)
304 {
305 char *buf;
306 laszip_get_error(m_p->laszip, &buf);
307 throwError(buf);
308 }
309 #endif
310 }
311
312
ready(PointTableRef table)313 void LasReader::ready(PointTableRef table)
314 {
315 createStream();
316 std::istream *stream(m_streamIf->m_istream);
317
318 m_p->index = 0;
319 if (m_p->header.compressed())
320 {
321 #ifdef PDAL_HAVE_LASZIP
322 if (m_args->compression == "LASZIP")
323 {
324 laszip_BOOL compressed;
325
326 handleLaszip(laszip_create(&m_p->laszip));
327 handleLaszip(laszip_open_reader_stream(m_p->laszip, *stream,
328 &compressed));
329 handleLaszip(laszip_get_point_pointer(m_p->laszip, &m_p->laszipPoint));
330 handleLaszip(laszip_seek_point(m_p->laszip, m_args->start));
331 }
332 #endif
333
334 #ifdef PDAL_HAVE_LAZPERF
335 if (m_args->compression == "LAZPERF")
336 {
337 delete m_p->decompressor;
338
339 const LasVLR *vlr = m_p->header.findVlr(LASZIP_USER_ID, LASZIP_RECORD_ID);
340 if (!vlr)
341 throwError("LAZ file missing required laszip VLR.");
342 int ebCount = m_p->header.pointLen() - m_p->header.basePointLen();
343 m_p->decompressor = new LazPerfVlrDecompressor(*stream, m_p->header.pointFormat(),
344 ebCount, m_p->header.pointOffset(), vlr->data());
345 if (m_args->start > 0)
346 {
347 if (m_args->start > m_p->header.pointCount())
348 throwError("'start' option set past end of file.");
349 m_p->decompressor->seek(m_args->start);
350 }
351 m_p->decompressorBuf.resize(m_p->header.pointLen());
352 }
353 #endif
354
355 #if !defined(PDAL_HAVE_LAZPERF) && !defined(PDAL_HAVE_LASZIP)
356 throwError("Can't read compressed file without LASzip or "
357 "LAZperf decompression library.");
358 #endif
359 }
360 else
361 {
362 std::istream::pos_type start = m_p->header.pointOffset() +
363 (m_args->start * m_p->header.pointLen());
364 stream->seekg(start);
365 }
366 }
367
368
369 namespace
370 {
371
addForwardMetadata(MetadataNode & forward,MetadataNode & m,const std::string & name,double val,const std::string description,size_t precision)372 void addForwardMetadata(MetadataNode& forward, MetadataNode& m,
373 const std::string& name, double val, const std::string description,
374 size_t precision)
375 {
376 MetadataNode n = m.add(name, val, description, precision);
377
378 // If the entry doesn't already exist, just add it.
379 MetadataNode f = forward.findChild(name);
380 if (!f.valid())
381 {
382 forward.add(n);
383 return;
384 }
385
386 // If the old value and new values aren't the same, set an invalid flag.
387 MetadataNode temp = f.addOrUpdate("temp", val, description, precision);
388 if (f.value<std::string>() != temp.value<std::string>())
389 forward.addOrUpdate(name + "INVALID", "");
390 }
391
392 }
393
394 // Store data in the normal metadata place. Also store it in the private
395 // lasforward metadata node.
396 template <typename T>
addForwardMetadata(MetadataNode & forward,MetadataNode & m,const std::string & name,T val,const std::string description)397 void addForwardMetadata(MetadataNode& forward, MetadataNode& m,
398 const std::string& name, T val, const std::string description)
399 {
400 MetadataNode n = m.add(name, val, description);
401
402 // If the entry doesn't already exist, just add it.
403 MetadataNode f = forward.findChild(name);
404 if (!f.valid())
405 {
406 forward.add(n);
407 return;
408 }
409
410 // If the old value and new values aren't the same, set an invalid flag.
411 MetadataNode temp = f.addOrUpdate("temp", val);
412 if (f.value<std::string>() != temp.value<std::string>())
413 forward.addOrUpdate(name + "INVALID", "");
414 }
415
416
extractHeaderMetadata(MetadataNode & forward,MetadataNode & m)417 void LasReader::extractHeaderMetadata(MetadataNode& forward, MetadataNode& m)
418 {
419 m.add<bool>("compressed", m_p->header.compressed(),
420 "true if this LAS file is compressed");
421
422 addForwardMetadata(forward, m, "major_version", m_p->header.versionMajor(),
423 "The major LAS version for the file, always 1 for now");
424 addForwardMetadata(forward, m, "minor_version", m_p->header.versionMinor(),
425 "The minor LAS version for the file");
426 addForwardMetadata(forward, m, "dataformat_id", m_p->header.pointFormat(),
427 "LAS Point Data Format");
428 if (m_p->header.versionAtLeast(1, 1))
429 addForwardMetadata(forward, m, "filesource_id",
430 m_p->header.fileSourceId(), "File Source ID (Flight Line Number "
431 "if this file was derived from an original flight line).");
432 if (m_p->header.versionAtLeast(1, 2))
433 {
434 // For some reason we've written global encoding as a base 64
435 // encoded value in the past. In an effort to standardize things,
436 // I'm writing this as a special value, and will also write
437 // global_encoding like we write all other header metadata.
438 uint16_t globalEncoding = m_p->header.globalEncoding();
439 m.addEncoded("global_encoding_base64", (uint8_t *)&globalEncoding,
440 sizeof(globalEncoding),
441 "Global Encoding: general property bit field.");
442
443 addForwardMetadata(forward, m, "global_encoding",
444 m_p->header.globalEncoding(),
445 "Global Encoding: general property bit field.");
446 }
447
448 addForwardMetadata(forward, m, "project_id", m_p->header.projectId(),
449 "Project ID.");
450 addForwardMetadata(forward, m, "system_id", m_p->header.systemId(),
451 "Generating system ID.");
452 addForwardMetadata(forward, m, "software_id", m_p->header.softwareId(),
453 "Generating software description.");
454 addForwardMetadata(forward, m, "creation_doy", m_p->header.creationDOY(),
455 "Day, expressed as an unsigned short, on which this file was created. "
456 "Day is computed as the Greenwich Mean Time (GMT) day. January 1 is "
457 "considered day 1.");
458 addForwardMetadata(forward, m, "creation_year", m_p->header.creationYear(),
459 "The year, expressed as a four digit number, in which the file was "
460 "created.");
461 addForwardMetadata(forward, m, "scale_x", m_p->header.scaleX(),
462 "The scale factor for X values.", 15);
463 addForwardMetadata(forward, m, "scale_y", m_p->header.scaleY(),
464 "The scale factor for Y values.", 15);
465 addForwardMetadata(forward, m, "scale_z", m_p->header.scaleZ(),
466 "The scale factor for Z values.", 15);
467 addForwardMetadata(forward, m, "offset_x", m_p->header.offsetX(),
468 "The offset for X values.", 15);
469 addForwardMetadata(forward, m, "offset_y", m_p->header.offsetY(),
470 "The offset for Y values.", 15);
471 addForwardMetadata(forward, m, "offset_z", m_p->header.offsetZ(),
472 "The offset for Z values.", 15);
473
474 m.add("point_length", m_p->header.pointLen(),
475 "The size, in bytes, of each point records.");
476 m.add("header_size", m_p->header.vlrOffset(),
477 "The size, in bytes, of the header block, including any extension "
478 "by specific software.");
479 m.add("dataoffset", m_p->header.pointOffset(),
480 "The actual number of bytes from the beginning of the file to the "
481 "first field of the first point record data field. This data offset "
482 "must be updated if any software adds data from the Public Header "
483 "Block or adds/removes data to/from the Variable Length Records.");
484 m.add<double>("minx", m_p->header.minX(),
485 "The max and min data fields are the actual unscaled extents of the "
486 "LAS point file data, specified in the coordinate system of the LAS "
487 "data.");
488 m.add<double>("miny", m_p->header.minY(),
489 "The max and min data fields are the actual unscaled extents of the "
490 "LAS point file data, specified in the coordinate system of the LAS "
491 "data.");
492 m.add<double>("minz", m_p->header.minZ(),
493 "The max and min data fields are the actual unscaled extents of the "
494 "LAS point file data, specified in the coordinate system of the LAS "
495 "data.");
496 m.add<double>("maxx", m_p->header.maxX(),
497 "The max and min data fields are the actual unscaled extents of the "
498 "LAS point file data, specified in the coordinate system of the LAS "
499 "data.");
500 m.add<double>("maxy", m_p->header.maxY(),
501 "The max and min data fields are the actual unscaled extents of the "
502 "LAS point file data, specified in the coordinate system of the LAS "
503 "data.");
504 m.add<double>("maxz", m_p->header.maxZ(),
505 "The max and min data fields are the actual unscaled extents of the "
506 "LAS point file data, specified in the coordinate system of the LAS "
507 "data.");
508 m.add<point_count_t>("count",
509 m_p->header.pointCount(), "This field contains the total "
510 "number of point records within the file.");
511
512 m.add<std::string>("gtiff", m_p->header.geotiffPrint(),
513 "GTifPrint output of GeoTIFF keys");
514
515 // PDAL metadata VLR
516 const LasVLR *vlr = m_p->header.findVlr("PDAL", 12);
517 if (vlr)
518 {
519 const char *pos = vlr->data();
520 size_t size = vlr->dataLen();
521 m.addWithType("pdal_metadata", std::string(pos, size), "json",
522 "PDAL Processing Metadata");
523 }
524 //
525 // PDAL pipeline VLR
526 vlr = m_p->header.findVlr("PDAL", 13);
527 if (vlr)
528 {
529 const char *pos = vlr->data();
530 size_t size = vlr->dataLen();
531 m.addWithType("pdal_pipeline", std::string(pos, size), "json",
532 "PDAL Processing Pipeline");
533 }
534 }
535
536
readExtraBytesVlr()537 void LasReader::readExtraBytesVlr()
538 {
539 const LasVLR *vlr = m_p->header.findVlr(SPEC_USER_ID,
540 EXTRA_BYTES_RECORD_ID);
541 if (!vlr)
542 return;
543 const char *pos = vlr->data();
544 size_t size = vlr->dataLen();
545 if (size % sizeof(ExtraBytesSpec) != 0)
546 {
547 log()->get(LogLevel::Warning) << "Bad size for extra bytes VLR. "
548 "Ignoring.";
549 return;
550 }
551 size /= sizeof(ExtraBytesSpec);
552 std::vector<ExtraBytesIf> ebList;
553 while (size--)
554 {
555 ExtraBytesIf eb;
556 eb.readFrom(pos);
557 ebList.push_back(eb);
558 pos += sizeof(ExtraBytesSpec);
559 }
560
561 std::vector<ExtraDim> extraDims;
562 for (ExtraBytesIf& eb : ebList)
563 {
564 std::vector<ExtraDim> eds = eb.toExtraDims();
565 for (auto& ed : eds)
566 extraDims.push_back(std::move(ed));
567 }
568 if (m_p->extraDims.size() && m_p->extraDims != extraDims)
569 log()->get(LogLevel::Warning) << "Extra byte dimensions specified "
570 "in pipeline and VLR don't match. Ignoring pipeline-specified "
571 "dimensions";
572 m_p->extraDims = extraDims;
573 }
574
575
setSrs(MetadataNode & m)576 void LasReader::setSrs(MetadataNode& m)
577 {
578 setSpatialReference(m, m_p->header.srs());
579 }
580
581
extractVlrMetadata(MetadataNode & forward,MetadataNode & m)582 void LasReader::extractVlrMetadata(MetadataNode& forward, MetadataNode& m)
583 {
584 static const size_t DATA_LEN_MAX = 1000000;
585
586 int i = 0;
587 for (auto vlr : m_p->header.vlrs())
588 {
589 if (vlr.dataLen() > DATA_LEN_MAX)
590 continue;
591
592 std::ostringstream name;
593 name << "vlr_" << i++;
594 MetadataNode vlrNode(name.str());
595
596 vlrNode.addEncoded("data",
597 (const uint8_t *)vlr.data(), vlr.dataLen(), vlr.description());
598 vlrNode.add("user_id", vlr.userId(),
599 "User ID of the record or pre-defined value from the "
600 "specification.");
601 vlrNode.add("record_id", vlr.recordId(),
602 "Record ID specified by the user.");
603 vlrNode.add("description", vlr.description());
604 m.add(vlrNode);
605
606 if (vlr.userId() == TRANSFORM_USER_ID||
607 vlr.userId() == LASZIP_USER_ID ||
608 vlr.userId() == LIBLAS_USER_ID)
609 continue;
610 if (vlr.userId() == SPEC_USER_ID &&
611 vlr.recordId() != 0 && vlr.recordId() != 3)
612 continue;
613 forward.add(vlrNode);
614 }
615 }
616
617
addDimensions(PointLayoutPtr layout)618 void LasReader::addDimensions(PointLayoutPtr layout)
619 {
620 using namespace Dimension;
621
622 layout->registerDim(Id::X, Type::Double);
623 layout->registerDim(Id::Y, Type::Double);
624 layout->registerDim(Id::Z, Type::Double);
625 layout->registerDim(Id::Intensity, Type::Unsigned16);
626 layout->registerDim(Id::ReturnNumber, Type::Unsigned8);
627 layout->registerDim(Id::NumberOfReturns, Type::Unsigned8);
628 layout->registerDim(Id::ScanDirectionFlag, Type::Unsigned8);
629 layout->registerDim(Id::EdgeOfFlightLine, Type::Unsigned8);
630 layout->registerDim(Id::Classification, Type::Unsigned8);
631 layout->registerDim(Id::ScanAngleRank, Type::Float);
632 layout->registerDim(Id::UserData, Type::Unsigned8);
633 layout->registerDim(Id::PointSourceId, Type::Unsigned16);
634
635 if (m_p->header.hasTime())
636 layout->registerDim(Id::GpsTime, Type::Double);
637 if (m_p->header.hasColor())
638 {
639 layout->registerDim(Id::Red, Type::Unsigned16);
640 layout->registerDim(Id::Green, Type::Unsigned16);
641 layout->registerDim(Id::Blue, Type::Unsigned16);
642 }
643 if (m_p->header.hasInfrared())
644 layout->registerDim(Id::Infrared);
645 if (m_p->header.has14PointFormat())
646 {
647 layout->registerDim(Id::ScanChannel);
648 layout->registerDim(Id::ClassFlags);
649 }
650
651 size_t ebLen = m_p->header.pointLen() - m_p->header.basePointLen();
652 for (auto& dim : m_p->extraDims)
653 {
654 if (dim.m_size > ebLen)
655 throwError("Extra byte specification exceeds point length beyond base format length.");
656 ebLen -= dim.m_size;
657
658 Dimension::Type type = dim.m_dimType.m_type;
659 if (type == Dimension::Type::None)
660 continue;
661 if (dim.m_dimType.m_xform.nonstandard())
662 type = Dimension::Type::Double;
663 if (m_args->fixNames)
664 dim.m_name = Dimension::fixName(dim.m_name);
665 dim.m_dimType.m_id = layout->registerOrAssignDim(dim.m_name, type);
666 }
667 }
668
669
processOne(PointRef & point)670 bool LasReader::processOne(PointRef& point)
671 {
672 if (m_p->index >= getNumPoints())
673 return false;
674
675 size_t pointLen = m_p->header.pointLen();
676
677 if (m_p->header.compressed())
678 {
679 #ifdef PDAL_HAVE_LASZIP
680 if (m_args->compression == "LASZIP")
681 {
682 handleLaszip(laszip_read_point(m_p->laszip));
683 loadPoint(point);
684 }
685 #endif
686
687 #ifdef PDAL_HAVE_LAZPERF
688 if (m_args->compression == "LAZPERF")
689 {
690 m_p->decompressor->decompress(m_p->decompressorBuf.data());
691 loadPoint(point, m_p->decompressorBuf.data(), pointLen);
692 }
693 #endif
694 #if !defined(PDAL_HAVE_LAZPERF) && !defined(PDAL_HAVE_LASZIP)
695 throwError("Can't read compressed file without LASzip or "
696 "LAZperf decompression library.");
697 #endif
698 } // compression
699 else
700 {
701 std::vector<char> buf(m_p->header.pointLen());
702
703 m_streamIf->m_istream->read(buf.data(), pointLen);
704 loadPoint(point, buf.data(), pointLen);
705 }
706 m_p->index++;
707 return true;
708 }
709
710
read(PointViewPtr view,point_count_t count)711 point_count_t LasReader::read(PointViewPtr view, point_count_t count)
712 {
713 size_t pointLen = m_p->header.pointLen();
714 count = (std::min)(count, getNumPoints() - m_p->index);
715
716 PointId i = 0;
717 if (m_p->header.compressed())
718 {
719 #if defined(PDAL_HAVE_LAZPERF) || defined(PDAL_HAVE_LASZIP)
720 if (m_args->compression == "LASZIP" || m_args->compression == "LAZPERF")
721 {
722 for (i = 0; i < count; i++)
723 {
724 PointRef point = view->point(i);
725 PointId id = view->size();
726 processOne(point);
727 if (m_cb)
728 m_cb(*view, id);
729 }
730 }
731 #else
732 throwError("Can't read compressed file without LASzip or "
733 "LAZperf decompression library.");
734 #endif
735 }
736 else
737 {
738 point_count_t remaining = count;
739
740 // Make a buffer at most a meg.
741 size_t bufsize = (std::min)((point_count_t)1000000, count * pointLen);
742 std::vector<char> buf(bufsize);
743 try
744 {
745 do
746 {
747 point_count_t blockPoints = readFileBlock(buf, remaining);
748 remaining -= blockPoints;
749 char *pos = buf.data();
750 while (blockPoints--)
751 {
752 PointId id = view->size();
753 PointRef point = view->point(id);
754 loadPoint(point, pos, pointLen);
755 if (m_cb)
756 m_cb(*view, id);
757 pos += pointLen;
758 i++;
759 }
760 } while (remaining);
761 }
762 catch (std::out_of_range&)
763 {}
764 catch (invalid_stream&)
765 {}
766 }
767 m_p->index += i;
768 return (point_count_t)i;
769 }
770
771
readFileBlock(std::vector<char> & buf,point_count_t maxpoints)772 point_count_t LasReader::readFileBlock(std::vector<char>& buf,
773 point_count_t maxpoints)
774 {
775 std::istream *stream(m_streamIf->m_istream);
776
777 size_t ptLen = m_p->header.pointLen();
778 point_count_t blockpoints = buf.size() / ptLen;
779
780 blockpoints = (std::min)(maxpoints, blockpoints);
781 if (stream->eof())
782 throw invalid_stream("stream is done");
783
784 stream->read(buf.data(), blockpoints * ptLen);
785 if (stream->gcount() != (std::streamsize)(blockpoints * ptLen))
786 {
787 // we read fewer bytes than we asked for
788 // because the file was either truncated
789 // or the header is bunk.
790 blockpoints = stream->gcount() / ptLen;
791 }
792 return blockpoints;
793 }
794
795
796 #ifdef PDAL_HAVE_LASZIP
loadPoint(PointRef & point)797 void LasReader::loadPoint(PointRef& point)
798 {
799 if (m_p->header.has14PointFormat())
800 loadPointV14(point);
801 else
802 loadPointV10(point);
803 }
804 #endif // PDAL_HAVE_LASZIP
805
806
loadPoint(PointRef & point,char * buf,size_t bufsize)807 void LasReader::loadPoint(PointRef& point, char *buf, size_t bufsize)
808 {
809 if (m_p->header.has14PointFormat())
810 loadPointV14(point, buf, bufsize);
811 else
812 loadPointV10(point, buf, bufsize);
813 }
814
815
816 #ifdef PDAL_HAVE_LASZIP
loadPointV10(PointRef & point)817 void LasReader::loadPointV10(PointRef& point)
818 {
819 // We used to pass the laszip point as an argument, but this allows us to keep
820 // any laszip information out of LasReader.hpp.
821 laszip_point& p = *m_p->laszipPoint;
822 const LasHeader& h = m_p->header;
823
824 double x = p.X * h.scaleX() + h.offsetX();
825 double y = p.Y * h.scaleY() + h.offsetY();
826 double z = p.Z * h.scaleZ() + h.offsetZ();
827
828 point.setField(Dimension::Id::X, x);
829 point.setField(Dimension::Id::Y, y);
830 point.setField(Dimension::Id::Z, z);
831 point.setField(Dimension::Id::Intensity, p.intensity);
832 point.setField(Dimension::Id::ReturnNumber, p.return_number);
833 point.setField(Dimension::Id::NumberOfReturns, p.number_of_returns);
834 point.setField(Dimension::Id::ScanDirectionFlag, p.scan_direction_flag);
835 point.setField(Dimension::Id::EdgeOfFlightLine, p.edge_of_flight_line);
836 uint8_t classification = p.classification | (p.synthetic_flag << 5) |
837 (p.keypoint_flag << 6) | (p.withheld_flag << 7);
838 point.setField(Dimension::Id::Classification, classification);
839 point.setField(Dimension::Id::ScanAngleRank, p.scan_angle_rank);
840 point.setField(Dimension::Id::UserData, p.user_data);
841 point.setField(Dimension::Id::PointSourceId, p.point_source_ID);
842
843 if (h.hasTime())
844 point.setField(Dimension::Id::GpsTime, p.gps_time);
845
846 if (h.hasColor())
847 {
848 point.setField(Dimension::Id::Red, p.rgb[0]);
849 point.setField(Dimension::Id::Green, p.rgb[1]);
850 point.setField(Dimension::Id::Blue, p.rgb[2]);
851 }
852
853 if (m_p->extraDims.size())
854 {
855 LeExtractor extractor((const char *)p.extra_bytes, p.num_extra_bytes);
856 loadExtraDims(extractor, point);
857 }
858 }
859 #endif // PDAL_HAVE_LASZIP
860
loadPointV10(PointRef & point,char * buf,size_t bufsize)861 void LasReader::loadPointV10(PointRef& point, char *buf, size_t bufsize)
862 {
863 LeExtractor istream(buf, bufsize);
864
865 int32_t xi, yi, zi;
866 istream >> xi >> yi >> zi;
867
868 const LasHeader& h = m_p->header;
869
870 double x = xi * h.scaleX() + h.offsetX();
871 double y = yi * h.scaleY() + h.offsetY();
872 double z = zi * h.scaleZ() + h.offsetZ();
873
874 uint16_t intensity;
875 uint8_t flags;
876 uint8_t classification;
877 int8_t scanAngleRank;
878 uint8_t user;
879 uint16_t pointSourceId;
880
881 istream >> intensity >> flags >> classification >> scanAngleRank >>
882 user >> pointSourceId;
883
884 uint8_t returnNum = flags & 0x07;
885 uint8_t numReturns = (flags >> 3) & 0x07;
886 uint8_t scanDirFlag = (flags >> 6) & 0x01;
887 uint8_t flight = (flags >> 7) & 0x01;
888
889 point.setField(Dimension::Id::X, x);
890 point.setField(Dimension::Id::Y, y);
891 point.setField(Dimension::Id::Z, z);
892 point.setField(Dimension::Id::Intensity, intensity);
893 point.setField(Dimension::Id::ReturnNumber, returnNum);
894 point.setField(Dimension::Id::NumberOfReturns, numReturns);
895 point.setField(Dimension::Id::ScanDirectionFlag, scanDirFlag);
896 point.setField(Dimension::Id::EdgeOfFlightLine, flight);
897 point.setField(Dimension::Id::Classification, classification);
898 point.setField(Dimension::Id::ScanAngleRank, scanAngleRank);
899 point.setField(Dimension::Id::UserData, user);
900 point.setField(Dimension::Id::PointSourceId, pointSourceId);
901
902 if (h.hasTime())
903 {
904 double time;
905 istream >> time;
906 point.setField(Dimension::Id::GpsTime, time);
907 }
908
909 if (h.hasColor())
910 {
911 uint16_t red, green, blue;
912 istream >> red >> green >> blue;
913 point.setField(Dimension::Id::Red, red);
914 point.setField(Dimension::Id::Green, green);
915 point.setField(Dimension::Id::Blue, blue);
916 }
917
918 if (m_p->extraDims.size())
919 loadExtraDims(istream, point);
920 }
921
922
923 #ifdef PDAL_HAVE_LASZIP
loadPointV14(PointRef & point)924 void LasReader::loadPointV14(PointRef& point)
925 {
926 // We used to pass the laszip point as an argument, but this allows us to keep
927 // any laszip information out of LasReader.hpp.
928 laszip_point& p = *m_p->laszipPoint;
929 const LasHeader& h = m_p->header;
930
931 double x = p.X * h.scaleX() + h.offsetX();
932 double y = p.Y * h.scaleY() + h.offsetY();
933 double z = p.Z * h.scaleZ() + h.offsetZ();
934
935 point.setField(Dimension::Id::X, x);
936 point.setField(Dimension::Id::Y, y);
937 point.setField(Dimension::Id::Z, z);
938 point.setField(Dimension::Id::Intensity, p.intensity);
939 point.setField(Dimension::Id::ReturnNumber, p.extended_return_number);
940 point.setField(Dimension::Id::NumberOfReturns,
941 p.extended_number_of_returns);
942 point.setField(Dimension::Id::ClassFlags, p.extended_classification_flags);
943 point.setField(Dimension::Id::ScanChannel, p.extended_scanner_channel);
944 point.setField(Dimension::Id::ScanDirectionFlag, p.scan_direction_flag);
945 point.setField(Dimension::Id::EdgeOfFlightLine, p.edge_of_flight_line);
946 point.setField(Dimension::Id::Classification, p.extended_classification);
947 point.setField(Dimension::Id::ScanAngleRank, p.extended_scan_angle * .006);
948 point.setField(Dimension::Id::UserData, p.user_data);
949 point.setField(Dimension::Id::PointSourceId, p.point_source_ID);
950 point.setField(Dimension::Id::GpsTime, p.gps_time);
951
952 if (h.hasColor())
953 {
954 point.setField(Dimension::Id::Red, p.rgb[0]);
955 point.setField(Dimension::Id::Green, p.rgb[1]);
956 point.setField(Dimension::Id::Blue, p.rgb[2]);
957 }
958
959 if (h.hasInfrared())
960 {
961 point.setField(Dimension::Id::Infrared, p.rgb[3]);
962 }
963
964 if (m_p->extraDims.size())
965 {
966 LeExtractor extractor((const char *)p.extra_bytes, p.num_extra_bytes);
967 loadExtraDims(extractor, point);
968 }
969 }
970 #endif // PDAL_HAVE_LASZIP
971
972
loadPointV14(PointRef & point,char * buf,size_t bufsize)973 void LasReader::loadPointV14(PointRef& point, char *buf, size_t bufsize)
974 {
975 LeExtractor istream(buf, bufsize);
976
977 int32_t xi, yi, zi;
978 istream >> xi >> yi >> zi;
979
980 const LasHeader& h = m_p->header;
981
982 double x = xi * h.scaleX() + h.offsetX();
983 double y = yi * h.scaleY() + h.offsetY();
984 double z = zi * h.scaleZ() + h.offsetZ();
985
986 uint16_t intensity;
987 uint8_t returnInfo;
988 uint8_t flags;
989 uint8_t classification;
990 uint8_t user;
991 int16_t scanAngle;
992 uint16_t pointSourceId;
993 double gpsTime;
994
995 istream >> intensity >> returnInfo >> flags >> classification >> user >>
996 scanAngle >> pointSourceId >> gpsTime;
997
998 uint8_t returnNum = returnInfo & 0x0F;
999 uint8_t numReturns = (returnInfo >> 4) & 0x0F;
1000 uint8_t classFlags = flags & 0x0F;
1001 uint8_t scanChannel = (flags >> 4) & 0x03;
1002 uint8_t scanDirFlag = (flags >> 6) & 0x01;
1003 uint8_t flight = (flags >> 7) & 0x01;
1004
1005 point.setField(Dimension::Id::X, x);
1006 point.setField(Dimension::Id::Y, y);
1007 point.setField(Dimension::Id::Z, z);
1008 point.setField(Dimension::Id::Intensity, intensity);
1009 point.setField(Dimension::Id::ReturnNumber, returnNum);
1010 point.setField(Dimension::Id::NumberOfReturns, numReturns);
1011 point.setField(Dimension::Id::ClassFlags, classFlags);
1012 point.setField(Dimension::Id::ScanChannel, scanChannel);
1013 point.setField(Dimension::Id::ScanDirectionFlag, scanDirFlag);
1014 point.setField(Dimension::Id::EdgeOfFlightLine, flight);
1015 point.setField(Dimension::Id::Classification, classification);
1016 point.setField(Dimension::Id::ScanAngleRank, scanAngle * .006);
1017 point.setField(Dimension::Id::UserData, user);
1018 point.setField(Dimension::Id::PointSourceId, pointSourceId);
1019 point.setField(Dimension::Id::GpsTime, gpsTime);
1020
1021 if (h.hasColor())
1022 {
1023 uint16_t red, green, blue;
1024 istream >> red >> green >> blue;
1025 point.setField(Dimension::Id::Red, red);
1026 point.setField(Dimension::Id::Green, green);
1027 point.setField(Dimension::Id::Blue, blue);
1028 }
1029
1030 if (h.hasInfrared())
1031 {
1032 uint16_t nearInfraRed;
1033
1034 istream >> nearInfraRed;
1035 point.setField(Dimension::Id::Infrared, nearInfraRed);
1036 }
1037
1038 if (m_p->extraDims.size())
1039 loadExtraDims(istream, point);
1040 }
1041
1042
loadExtraDims(LeExtractor & istream,PointRef & point)1043 void LasReader::loadExtraDims(LeExtractor& istream, PointRef& point)
1044 {
1045 for (auto& dim : m_p->extraDims)
1046 {
1047 // Dimension type of None is undefined and unprocessed
1048 if (dim.m_dimType.m_type == Dimension::Type::None)
1049 {
1050 istream.skip(dim.m_size);
1051 continue;
1052 }
1053
1054 Everything e = Utils::extractDim(istream, dim.m_dimType.m_type);
1055 if (dim.m_dimType.m_xform.nonstandard())
1056 {
1057 double d = Utils::toDouble(e, dim.m_dimType.m_type);
1058 d = d * dim.m_dimType.m_xform.m_scale.m_val +
1059 dim.m_dimType.m_xform.m_offset.m_val;
1060 point.setField(dim.m_dimType.m_id, d);
1061 }
1062 else
1063 point.setField(dim.m_dimType.m_id, dim.m_dimType.m_type, &e);
1064 }
1065 }
1066
1067
done(PointTableRef)1068 void LasReader::done(PointTableRef)
1069 {
1070 #ifdef PDAL_HAVE_LASZIP
1071 if (m_p->laszip)
1072 {
1073 handleLaszip(laszip_close_reader(m_p->laszip));
1074 handleLaszip(laszip_destroy(m_p->laszip));
1075 }
1076 #endif
1077 m_streamIf.reset();
1078 }
1079
eof()1080 bool LasReader::eof()
1081 {
1082 return m_p->index >= getNumPoints();
1083 }
1084
1085
1086 } // namespace pdal
1087