Merge remote-tracking branch 'origin/master' into python-numpy-array-…

…direct

CMakeLists.txt

@@ -216,6 +216,7 @@ if (NOT PDAL_HAVE_LIBXML2)
     file(GLOB XML_SRCS
         io/Ilvis2MetadataReader.cpp
         io/Ilvis2Metadata.cpp
+        io/Ilvis2Reader.cpp
         ${PDAL_SRC_DIR}/DbWriter.cpp
         ${PDAL_SRC_DIR}/DbReader.cpp
         ${PDAL_SRC_DIR}/XMLSchema.cpp)

cmake/libxml2.cmake

@@ -12,6 +12,11 @@ mark_as_advanced(CLEAR LIBXML2_INCLUDE_DIR)
 mark_as_advanced(CLEAR LIBXML2_LIBRARIES)
 if (LIBXML2_FOUND)
     set(PDAL_HAVE_LIBXML2 1)
+else()
+    unset(LIBXML2_INCLUDE_DIR CACHE)
+    unset(LIBXML2_LIBRARY CACHE)
+    unset(LIBXML2_XMLLINT_EXECUTABLE CACHE)
+    unset(LIBXML2_LIBRARIES) # Find-module output variable, not cache variable
 endif()
 
 set_property(GLOBAL PROPERTY _LIBXML2_INCLUDED TRUE)

doc/stages/writers.pgpointcloud.rst

@@ -70,7 +70,6 @@ compression
 
   * **none** applies no compression
   * **dimensional** applies dynamic compression to each dimension separately
-  * **ght** applies a "geohash tree" compression by sorting the points into a prefix tree
   * **lazperf** applies a "laz" compression (using the `laz-perf`_ library in PostgreSQL Pointcloud)
 
 overwrite

doc/stages/writers.ply.rst

@@ -43,7 +43,7 @@ filename
 storage_mode
   Type of ply file to write. Valid values are 'ascii', 'little endian',
   'big endian', and 'default'. 'default' is binary output in the endianness
-  of the machine. [Default: 'default']
+  of the machine. [Default: 'ascii']
 
 dims
   List of dimensions to write as elements. [Default: all dimensions]
@@ -52,4 +52,11 @@ faces
   Write a mesh as faces in addition to writing points as vertices.
   [Default: false]
 
+precision
+  If specified, the number of digits to the right of the decimal place
+  using f-style formatting.  Only permitted when 'storage_mode' is 'ascii'.
+  See the `printf`_ reference for more information.
+  [Default: g-style formatting (variable precision)]
+
 .. _polygon file format: http://paulbourke.net/dataformats/ply/
+.. _printf: https://en.cppreference.com/w/cpp/io/c/fprintf

filters/ApproximateCoplanarFilter.cpp

@@ -79,8 +79,7 @@ void ApproximateCoplanarFilter::filter(PointView& view)
 {
     using namespace Eigen;
 
-    KD3Index kdi(view);
-    kdi.build();
+    KD3Index& kdi = view.build3dIndex();
 
     for (PointId i = 0; i < view.size(); ++i)
     {

filters/EigenvaluesFilter.cpp

@@ -78,8 +78,7 @@ void EigenvaluesFilter::filter(PointView& view)
 {
     using namespace Eigen;
 
-    KD3Index kdi(view);
-    kdi.build();
+    KD3Index& kdi = view.build3dIndex();
 
     for (PointId i = 0; i < view.size(); ++i)
     {

filters/EstimateRankFilter.cpp

@@ -73,8 +73,7 @@ void EstimateRankFilter::addDimensions(PointLayoutPtr layout)
 
 void EstimateRankFilter::filter(PointView& view)
 {
-    KD3Index kdi(view);
-    kdi.build();
+    KD3Index& kdi = view.build3dIndex();
 
     for (PointId i = 0; i < view.size(); ++i)
     {

filters/GreedyProjection.cpp

@@ -126,9 +126,9 @@ void GreedyProjection::addTriangle(PointId a, PointId b, PointId c)
 
 void GreedyProjection::filter(PointView& view)
 {
-    pdal::NormalFilter().doFilter(view);
+    NormalFilter().doFilter(view);
 
-    pdal::KD3Index& tree = view.build3dIndex();
+    KD3Index& tree = view.build3dIndex();
 
     view_ = &view;
     mesh_ = view_->createMesh(getName());
@@ -154,17 +154,17 @@ void GreedyProjection::filter(PointView& view)
     already_connected_ = false; // see declaration for comments :P
 
     // initializing states and fringe neighbors
-    part_.clear ();
-    state_.clear ();
-    source_.clear ();
-    ffn_.clear ();
-    sfn_.clear ();
-    part_.resize(view.size ()); // indices of point's part
-    state_.resize(view.size (), GP3Type::FREE);
-    source_.resize(view.size ());
+    part_.clear();
+    state_.clear();
+    source_.clear();
+    ffn_.clear();
+    sfn_.clear();
+    part_.resize(view.size()); // indices of point's part
+    state_.resize(view.size(), GP3Type::FREE);
+    source_.resize(view.size());
     ffn_.resize(view.size());
     sfn_.resize(view.size());
-    fringe_queue_.clear ();
+    fringe_queue_.clear();
     int fqIdx = 0; // current fringe's index in the queue to be processed
 
   // Avoiding NaN coordinates if needed

filters/HAGFilter.cpp

@@ -96,8 +96,7 @@ void HAGFilter::filter(PointView& view)
             "as ground");
 
     // Build the 2D KD-tree.
-    KD2Index kdi(*gView);
-    kdi.build();
+    KD2Index& kdi = gView->build2dIndex();
 
     // Second pass: Find Z difference between non-ground points and the nearest
     // neighbor (2D) in the ground view.

filters/KDistanceFilter.cpp

@@ -71,10 +71,7 @@ void KDistanceFilter::filter(PointView& view)
 {
     using namespace Dimension;
 
-    // Build the 3D KD-tree.
-    log()->get(LogLevel::Debug) << "Building 3D KD-tree...\n";
-    KD3Index index(view);
-    index.build();
+    KD3Index& index = view.build3dIndex();
 
     // Increment the minimum number of points, as knnSearch will be returning
     // the neighbors along with the query point.

filters/LOFFilter.cpp

@@ -73,10 +73,7 @@ void LOFFilter::filter(PointView& view)
 {
     using namespace Dimension;
 
-    // Build the 3D KD-tree.
-    KD3Index index(view);
-    log()->get(LogLevel::Debug) << "Building 3D KD-tree...\n";
-    index.build();
+    KD3Index& index = view.build3dIndex();
 
     // Increment the minimum number of points, as knnSearch will be returning
     // the neighbors along with the query point.

filters/NNDistanceFilter.cpp

@@ -108,7 +108,6 @@ void NNDistanceFilter::filter(PointView& view)
     using namespace Dimension;
 
     // Build the 3D KD-tree.
-    log()->get(LogLevel::Debug) << "Building 3D KD-tree...\n";
     KD3Index& index = view.build3dIndex();
 
     // Increment the minimum number of points, as knnSearch will be returning

filters/NeighborClassifierFilter.cpp

@@ -171,8 +171,7 @@ void NeighborClassifierFilter::filter(PointView& view)
     PointRef point_src(view, 0);
     if (m_candidateFile.empty())
     {   // No candidate file so NN comes from src file
-        KD3Index kdiSrc(view);
-        kdiSrc.build();
+        KD3Index& kdiSrc = view.build3dIndex();
         PointRef point_nn(view, 0);
         for (PointId id = 0; id < view.size(); ++id)
         {
@@ -184,8 +183,7 @@ void NeighborClassifierFilter::filter(PointView& view)
     {   // NN comes from candidate file
         PointTable candTable;
         PointViewPtr candView = loadSet(m_candidateFile, candTable);
-        KD3Index kdiCand(*candView);
-        kdiCand.build();
+        KD3Index& kdiCand = candView->build3dIndex();
         PointRef point_nn(*candView, 0);
         for (PointId id = 0; id < view.size(); ++id)
         {

filters/PMFFilter.cpp

@@ -230,8 +230,7 @@ void PMFFilter::processGround(PointViewPtr view)
     }
 
     // build the 2D KD-tree
-    KD2Index kdi(*temp);
-    kdi.build();
+    KD2Index& kdi = temp->build2dIndex();
 
     // loop through all cells, and for each NaN, replace with elevation of
     // nearest neighbor

filters/RadialDensityFilter.cpp

@@ -72,9 +72,7 @@ void RadialDensityFilter::filter(PointView& view)
     using namespace Dimension;
 
     // Build the 3D KD-tree.
-    log()->get(LogLevel::Debug) << "Building 3D KD-tree...\n";
-    KD3Index index(view);
-    index.build();
+    KD3Index& index = view.build3dIndex();
 
     // Search for neighboring points within the specified radius. The number of
     // neighbors (which includes the query point) is normalized by the volume

filters/SMRFilter.cpp

@@ -548,8 +548,7 @@ std::vector<double> SMRFilter::knnfill(PointViewPtr view,
         }
     }
 
-    KD2Index kdi(*temp);
-    kdi.build();
+    KD2Index& kdi = view->build2dIndex();
 
     // Where the raster has voids (i.e., NaN), we search for that cell's eight
     // nearest neighbors, and fill the void with the average value of the

filters/SampleFilter.cpp

@@ -83,8 +83,7 @@ PointViewSet SampleFilter::run(PointViewPtr inView)
     PointViewPtr outView = inView->makeNew();
 
     // Build the 3D KD-tree.
-    KD3Index index(*inView);
-    index.build();
+    KD3Index& index = inView->build3dIndex();
 
     // All points are marked as kept (1) by default. As they are masked by
     // neighbors within the user-specified radius, their value is changed to 0.

filters/VoxelCenterNearestNeighborFilter.cpp

@@ -68,36 +68,44 @@ PointViewSet VoxelCenterNearestNeighborFilter::run(PointViewPtr view)
     BOX3D bounds;
     view->calculateBounds(bounds);
 
-    KD3Index kdi(*view);
-    kdi.build();
-
-    // Make an initial pass through the input PointView to detect populated
-    // voxels.
-    std::set<std::tuple<size_t, size_t, size_t>> populated_voxels;
+    // Find distance from voxel center to point.  If the distance is less
+    // than previous (or is the first one for the voxel), store the
+    // point ID and distance.
+    std::map<std::tuple<size_t, size_t, size_t>,
+        std::tuple<PointId, double>> populated_voxels;
     for (PointId id = 0; id < view->size(); ++id)
     {
-        double y = view->getFieldAs<double>(Dimension::Id::Y, id);
         double x = view->getFieldAs<double>(Dimension::Id::X, id);
+        double y = view->getFieldAs<double>(Dimension::Id::Y, id);
         double z = view->getFieldAs<double>(Dimension::Id::Z, id);
-        size_t r = static_cast<size_t>((y - bounds.miny) / m_cell);
         size_t c = static_cast<size_t>((x - bounds.minx) / m_cell);
+        size_t r = static_cast<size_t>((y - bounds.miny) / m_cell);
         size_t d = static_cast<size_t>((z - bounds.minz) / m_cell);
-        populated_voxels.emplace(std::make_tuple(r, c, d));
+        double xv = bounds.minx + (c + 0.5) * m_cell;
+        double yv = bounds.miny + (r + 0.5) * m_cell;
+        double zv = bounds.minz + (d + 0.5) * m_cell;
+        double dist = pow(xv - x, 2) + pow(yv - y, 2) + pow(zv - z, 2);
+
+        auto t = std::make_tuple(r, c, d);
+        auto pi = populated_voxels.find(t);
+        if (pi == populated_voxels.end())
+            populated_voxels.insert(
+                std::make_pair(t, std::make_tuple(id, dist)));
+        else
+        {
+            auto& t2 = pi->second;  // Get point/distance tuple.
+            double curDist = std::get<1>(t2);
+            if (dist < curDist)
+                t2 = std::make_tuple(id, dist);
+        }
     }
 
-    // Make a second pass through the populated voxels to find the nearest
-    // neighbor to each voxel center.
+    // Append the ID of the point nearest the voxel center to the output view.
     PointViewPtr output = view->makeNew();
     for (auto const& t : populated_voxels)
     {
-        auto& r = std::get<0>(t);
-        auto& c = std::get<1>(t);
-        auto& d = std::get<2>(t);
-        double y = bounds.miny + (r + 0.5) * m_cell;
-        double x = bounds.minx + (c + 0.5) * m_cell;
-        double z = bounds.minz + (d + 0.5) * m_cell;
-        std::vector<PointId> neighbors = kdi.neighbors(x, y, z, 1);
-        output->appendPoint(*view, neighbors[0]);
+        auto& t2 = t.second; // Get point/distance tuple.
+        output->appendPoint(*view, std::get<0>(t2));
     }
 
     PointViewSet viewSet;

filters/VoxelCentroidNearestNeighborFilter.cpp

@@ -71,8 +71,7 @@ PointViewSet VoxelCentroidNearestNeighborFilter::run(PointViewPtr view)
     BOX3D bounds;
     view->calculateBounds(bounds);
 
-    KD3Index kdi(*view);
-    kdi.build();
+    KD3Index& kdi = view->build3dIndex();
 
     // Make an initial pass through the input PointView to index PointIds by
     // row, column, and depth.

io/PlyWriter.cpp

@@ -66,11 +66,15 @@ void PlyWriter::addArgs(ProgramArgs& args)
     args.add("storage_mode", "PLY Storage Mode", m_format, Format::Ascii);
     args.add("dims", "Dimension names", m_dimNames);
     args.add("faces", "Write faces", m_faces);
+    m_precisionArg = &args.add("precision", "Output precision", m_precision, 3);
 }
 
 
 void PlyWriter::prepared(PointTableRef table)
 {
+    if (m_precisionArg->set() && m_format != Format::Ascii)
+        throwError("Option 'precision' can only be set of the 'storage_mode' "
+            "is ascii.");
     if (m_dimNames.size())
     {
         for (auto& name : m_dimNames)
@@ -149,6 +153,11 @@ void PlyWriter::ready(PointTableRef table)
             " points supported.");
 
     m_stream = Utils::createFile(m_filename, true);
+    if (m_format == Format::Ascii && m_precisionArg->set())
+    {
+        *m_stream << std::fixed;
+        m_stream->precision(m_precision);
+    }
     writeHeader(table.layout());
 }
 

io/PlyWriter.hpp

@@ -73,6 +73,8 @@ class PDAL_DLL PlyWriter : public Writer
     bool m_faces;
     StringList m_dimNames;
     Dimension::IdList m_dims;
+    int m_precision;
+    Arg *m_precisionArg;
     std::vector<PointViewPtr> m_views;
 };
 

pdal/PointTable.cpp

@@ -125,6 +125,23 @@ char *PointTable::getPoint(PointId idx)
 }
 
 
+ContiguousPointTable::~ContiguousPointTable()
+{}
+
+
+PointId ContiguousPointTable::addPoint()
+{
+    m_buf.resize(pointsToBytes(m_numPts + 1));
+    return m_numPts++;
+}
+
+
+char *ContiguousPointTable::getPoint(PointId idx)
+{
+    return m_buf.data() + pointsToBytes(idx);
+}
+
+
 MetadataNode BasePointTable::toMetadata() const
 {
     return layout()->toMetadata();

pdal/PointTable.hpp

@@ -166,6 +166,28 @@ class PDAL_DLL PointTable : public SimplePointTable
     PointLayout m_layout;
 };
 
+class PDAL_DLL ContiguousPointTable : public SimplePointTable
+{
+private:
+    std::vector<char> m_buf;
+    point_count_t m_numPts;
+
+public:
+    ContiguousPointTable() : SimplePointTable(m_layout), m_numPts(0)
+        {}
+    virtual ~ContiguousPointTable();
+    virtual bool supportsView() const
+        { return true; }
+
+protected:
+    virtual char *getPoint(PointId idx);
+
+private:
+    virtual PointId addPoint();
+
+    PointLayout m_layout;
+};
+
 /// A StreamPointTable must provide storage for point data up to its capacity.
 /// It must implement getPoint() which returns a pointer to a buffer of
 /// sufficient size to contain a point's data.  The minimum size required

pdal/PointView.cpp

@@ -237,7 +237,7 @@ void PointView::dump(std::ostream& ostr) const
         {
             Dimension::Id d = *di;
             const Dimension::Detail *dd = layout->dimDetail(d);
-            ostr << Dimension::name(d) << " (" <<
+            ostr << layout->dimName(d) << " (" <<
                 Dimension::interpretationName(dd->type()) << ") : ";
 
             switch (dd->type())