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