From b88c0390fcd99a79e002faaa50e692c9a1f21049 Mon Sep 17 00:00:00 2001 From: chambbj Date: Tue, 14 May 2019 15:26:26 -0500 Subject: [PATCH 1/6] Fix the broken Conda links (#2514) --- doc/quickstart.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/quickstart.rst b/doc/quickstart.rst index df191c236b..710579304e 100644 --- a/doc/quickstart.rst +++ b/doc/quickstart.rst @@ -38,9 +38,9 @@ Install Conda Conda installation instructions can be found at the following links. Read through them a bit for your platform so you have an idea what to expect. -* `Windows `__ -* `macOS `__ -* `Linux `__ +* `Windows `__ +* `macOS `__ +* `Linux `__ .. note:: @@ -56,7 +56,7 @@ On macOS and Linux, all Conda commands are typed into a terminal window. On Windows, commands are typed into the Anaconda Prompt window. Instructions can be found in the Conda `Getting Started`_ guide. -.. _`Getting Started`: https://conda.io/docs/user-guide/getting-started.html#starting-conda +.. _`Getting Started`: https://conda.io/projects/conda/en/latest/user-guide/getting-started.html#starting-conda Test Installation @@ -79,7 +79,7 @@ the PDAL maintenance branch. It is a good idea to install PDAL in it's own environment (or add it to an existing one). You will **NOT** want to add it to your default environment named ``base``. Managing environments is beyond the scope of - the quickstart, but you can read more about it `here `_. + the quickstart, but you can read more about it `here `_. To install the PDAL package so that we can use it to run PDAL commands, we run the following command to create an environment named ``myenv``, installing PDAL From f99564b6fb9a4b66b47cd8f301f39bbeb534d92c Mon Sep 17 00:00:00 2001 From: chambbj Date: Wed, 15 May 2019 10:58:39 -0500 Subject: [PATCH 2/6] Remove Mongus filter and all references (#2516) * Remove Mongus filter and all references * One more Mongus ref to be removed --- doc/references.rst | 2 - doc/stages/filters.mongus.rst | 64 ---- doc/stages/filters.rst | 4 - filters/MongusFilter.cpp | 567 ---------------------------------- filters/MongusFilter.hpp | 89 ------ 5 files changed, 726 deletions(-) delete mode 100644 doc/stages/filters.mongus.rst delete mode 100644 filters/MongusFilter.cpp delete mode 100644 filters/MongusFilter.hpp diff --git a/doc/references.rst b/doc/references.rst index 43c2249f23..70242c2ab6 100644 --- a/doc/references.rst +++ b/doc/references.rst @@ -52,8 +52,6 @@ Reference .. [Mesh2009] ALoopingIcon. "Meshing Point Clouds." *MESHLAB STUFF*. n.p., 7 Sept. 2009. Web. 13 Nov. 2015. -.. [Mongus2012] Mongus, D., Zalik, B., 2012. Parameter-free ground filtering of LiDAR data for automatic DTM generation. ISPRS J. Photogramm. Remote Sens. 67, 1–12. - .. [Pingel2013] Pingel, Thomas J., Keith C. Clarke, and William A. McBride. β€œAn Improved Simple Morphological Filter for the Terrain Classification of Airborne LIDAR Data.” ISPRS Journal of Photogrammetry and Remote Sensing 77 (2013): 21–30. .. [Rusu2008] Rusu, Radu Bogdan, et al. "Towards 3D point cloud based object maps for household environments." Robotics and Autonomous Systems 56.11 (2008): 927-941. diff --git a/doc/stages/filters.mongus.rst b/doc/stages/filters.mongus.rst deleted file mode 100644 index 83236dc7b5..0000000000 --- a/doc/stages/filters.mongus.rst +++ /dev/null @@ -1,64 +0,0 @@ -.. _filters.mongus: - -filters.mongus -=============================================================================== - -The **mongus filter** determines ground returns using the approach -outlined in [Mongus2012]_. - -.. note:: - - The mongus filter is deprecated and has been replaced by - :ref:`filters.pmf` and :ref:`filters.smrf`. - -The current implementation of this filter differs slightly from the -original paper. We weren't too happy with the criteria for how control -points at -the current level are compared against the TPS at the previous scale and were -exploring some alternate metrics. - -Some warts about the current implementation: - -* It writes many intermediate/debugging outputs to the current directory - while processing. - -* We require the specification of a max level, whereas the original paper - automatically determined an appropriate max level. - -.. embed:: - -Example -------- - -The sample pipeline below uses the filter to segment ground and -non-ground returns, writing only the ground returns to the output file. - -.. code-block:: json - - [ - "input.las", - { - "type":"filters.mongus", - "extract":true - }, - "output.laz" - ] - -Options -------------------------------------------------------------------------------- - -cell - Cell size. [Default: 1.0] - -classify - Apply classification labels (i.e., ground = 2)? [Default: true] - -extract - Extract ground returns (non-ground returns are cropped)? [Default: false] - -k - Standard deviation multiplier to be used when thresholding - values. [Default: 3.0] - -l - Maximum level in the hierarchical decomposition. [Default: 8] diff --git a/doc/stages/filters.rst b/doc/stages/filters.rst index 8a131ef5a3..4133f10102 100644 --- a/doc/stages/filters.rst +++ b/doc/stages/filters.rst @@ -45,7 +45,6 @@ invalidate an existing KD-tree. filters.info filters.kdistance filters.lof - filters.mongus filters.neighborclassifier filters.nndistance filters.normal @@ -98,9 +97,6 @@ invalidate an existing KD-tree. Compute pointwise Local Outlier Factor (along with K-Distance and Local Reachability Distance). -:ref:`filters.mongus` - Label ground/non-ground returns using [Mongus2012]_. - :ref:`filters.neighborclassifier` Update pointwise classification using k-nearest neighbor consensus voting. diff --git a/filters/MongusFilter.cpp b/filters/MongusFilter.cpp deleted file mode 100644 index 7ee628f473..0000000000 --- a/filters/MongusFilter.cpp +++ /dev/null @@ -1,567 +0,0 @@ -/****************************************************************************** -* Copyright (c) 2016, Bradley J Chambers (brad.chambers@gmail.com) -* -* All rights reserved. -* -* Redistribution and use in source and binary forms, with or without -* modification, are permitted provided that the following -* conditions are met: -* -* * Redistributions of source code must retain the above copyright -* notice, this list of conditions and the following disclaimer. -* * Redistributions in binary form must reproduce the above copyright -* notice, this list of conditions and the following disclaimer in -* the documentation and/or other materials provided -* with the distribution. -* * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the -* names of its contributors may be used to endorse or promote -* products derived from this software without specific prior -* written permission. -* -* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, -* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, -* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS -* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED -* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT -* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY -* OF SUCH DAMAGE. -****************************************************************************/ - -#include "MongusFilter.hpp" - -#include -#include -#include -#include -#include -#include - -#include - -// NOTE: This filter is deprecated as of 3/2019 - -namespace pdal -{ - -static StaticPluginInfo const s_info -{ - "filters.mongus", - "Mongus and Zalik (2012)", - "http://pdal.io/stages/filters.mongus.html" -}; - -CREATE_STATIC_STAGE(MongusFilter, s_info) - -std::string MongusFilter::getName() const -{ - return s_info.name; -} - -void MongusFilter::addArgs(ProgramArgs& args) -{ - args.add("cell", "Cell size", m_cellSize, 1.0); - args.add("k", "Stdev multiplier for threshold", m_k, 3.0); - args.add("l", "Max level", m_l, 8); - args.add("classify", "Apply classification labels?", m_classify, true); - args.add("extract", "Extract ground returns?", m_extract); -} - -void MongusFilter::addDimensions(PointLayoutPtr layout) -{ - layout->registerDim(Dimension::Id::Classification); -} - -int MongusFilter::getColIndex(double x, double cell_size) -{ - return static_cast(floor((x - m_bounds.minx) / cell_size)); -} - -int MongusFilter::getRowIndex(double y, double cell_size) -{ - return static_cast(floor((m_maxRow - y) / cell_size)); -} - -void MongusFilter::writeControl(Eigen::MatrixXd cx, Eigen::MatrixXd cy, Eigen::MatrixXd cz, std::string filename) -{ - using namespace Dimension; - - PipelineManager m; - - PointTable table; - PointViewPtr view(new PointView(table)); - - table.layout()->registerDim(Id::X); - table.layout()->registerDim(Id::Y); - table.layout()->registerDim(Id::Z); - - PointId i = 0; - for (auto j = 0; j < cz.size(); ++j) - { - if (std::isnan(cx(j)) || std::isnan(cy(j)) || std::isnan(cz(j))) - continue; - view->setField(Id::X, i, cx(j)); - view->setField(Id::Y, i, cy(j)); - view->setField(Id::Z, i, cz(j)); - i++; - } - - BufferReader r; - r.addView(view); - - Stage& w = m.makeWriter(filename, "writers.las", r); - w.prepare(table); - w.execute(table); -} - -std::vector MongusFilter::processGround(PointViewPtr view) -{ - using namespace Eigen; - - point_count_t np(view->size()); - - std::vector groundIdx; - - // initialization - - view->calculateBounds(m_bounds); - SpatialReference srs(view->spatialReference()); - - m_numCols = - static_cast(ceil((m_bounds.maxx - m_bounds.minx)/m_cellSize)) + 1; - m_numRows = - static_cast(ceil((m_bounds.maxy - m_bounds.miny)/m_cellSize)) + 1; - m_maxRow = static_cast(m_bounds.miny + m_numRows * m_cellSize); - - // create control points matrix at default cell size - MatrixXd cx(m_numRows, m_numCols); - cx.setConstant(std::numeric_limits::quiet_NaN()); - - MatrixXd cy(m_numRows, m_numCols); - cy.setConstant(std::numeric_limits::quiet_NaN()); - - MatrixXd cz(m_numRows, m_numCols); - cz.setConstant((std::numeric_limits::max)()); - - // find initial set of Z minimums at native resolution - for (point_count_t i = 0; i < np; ++i) - { - using namespace Dimension; - double x = view->getFieldAs(Id::X, i); - double y = view->getFieldAs(Id::Y, i); - double z = view->getFieldAs(Id::Z, i); - - int c = Utils::clamp(getColIndex(x, m_cellSize), 0, m_numCols-1); - int r = Utils::clamp(getRowIndex(y, m_cellSize), 0, m_numRows-1); - - if (z < cz(r, c)) - { - cx(r, c) = x; - cy(r, c) = y; - cz(r, c) = z; - } - } - - writeControl(cx, cy, cz, "grid_mins.laz"); - - // In our case, 2D structural elements of circular shape are employed and - // sufficient accuracy is achieved by using a larger window size for opening - // (W11) than for closing (W9). - MatrixXd mo = eigen::matrixOpen(cz, 11); - writeControl(cx, cy, mo, "grid_open.laz"); - MatrixXd mc = eigen::matrixClose(mo, 9); - writeControl(cx, cy, mc, "grid_close.laz"); - - // ...in order to minimize the distortions caused by such filtering, the - // output points ... are compared to C and only ci with significantly lower - // elevation [are] replaced... In our case, d = 1.0 m was used. - for (auto i = 0; i < cz.size(); ++i) - { - if ((mc(i) - cz(i)) >= 1.0) - cz(i) = mc(i); - } - // cz is still at native resolution, with low points replaced by - // morphological operators - writeControl(cx, cy, cz, "grid_mins_adjusted.laz"); - - // downsample control at max_level - int level = m_l; - double cur_cell_size = m_cellSize * std::pow(2, level); - // for max level = 8 and cell size 1, this is 256 - - MatrixXd x_prev, y_prev, z_prev; - - // Top-level control samples are assumed to be ground points, no filtering - // is applied. - downsampleMin(&cx, &cy, &cz, &x_prev, &y_prev, &z_prev, cur_cell_size); - // x|y|z_prev are control points downsampled to coarsest resolution for - // the hierarchy, e.g., for 512x512, this would be 2x2 - writeControl(x_prev, y_prev, z_prev, "control_init.laz"); - - // Point-filtering is performed iteratively at each level of the - // control-points hierarchy in a top-down fashion - for (auto l = level-1; l > 0; --l) - { - std::cerr << "Level " << l << std::endl; - cur_cell_size /= 2; - // 128, 64, 32, 16, 8, 4, 1 - - // compute TPS with update control at level - - // The interpolated surface is estimated based on the filtered set of - // TPS control-points at the previous level of hierarchy - // MatrixXd surface = TPS(x_prev, y_prev, z_prev, cur_cell_size); - // 4x4, 8x8, 16x16, 32x32, 64x64, 128x128, 256x256 - - // downsample control at level - MatrixXd x_samp, y_samp, z_samp; - downsampleMin(&cx, &cy, &cz, &x_samp, &y_samp, &z_samp, cur_cell_size); - // 4x4, 8x8, 16x16, 32x32, 64x64, 128x128, 256x256 - - MatrixXd surface = eigen::computeSpline(x_prev, y_prev, z_prev, - x_samp, y_samp); - - // if (l == 3) - // { - // log()->get(LogLevel::Debug) << cx.rows() << "\t" << cx.cols() << std::endl; - // log()->get(LogLevel::Debug) << x_prev.rows() << "\t" << x_prev.cols() << std::endl; - // log()->get(LogLevel::Debug) << x_samp.rows() << "\t" << x_samp.cols() << std::endl; - // log()->get(LogLevel::Debug) << surface.rows() << "\t" << surface.cols() << std::endl; - // log()->get(LogLevel::Debug) << "x: " << cx.row(1) << std::endl; - // log()->get(LogLevel::Debug) << "z: " << cz.row(1) << std::endl; - // log()->get(LogLevel::Debug) << "control_x: " << x_prev.row(0) << std::endl; - // log()->get(LogLevel::Debug) << "control_z: " << z_prev.row(0) << std::endl; - // log()->get(LogLevel::Debug) << "samples_x: " << x_samp.row(0) << std::endl; - // log()->get(LogLevel::Debug) << "samples_z: " << z_samp.row(0) << std::endl; - // log()->get(LogLevel::Debug) << "spline: " << surface.row(0) << std::endl; - // } - - char bufs[256]; - sprintf(bufs, "cur_control_%d.laz", l); - std::string names(bufs); - writeControl(x_samp, y_samp, z_samp, names); - - MatrixXd R = z_samp - surface; - - if (l == 7) - log()->get(LogLevel::Debug) << R << std::endl; - - double sum = 0.0; - double maxcoeff = std::numeric_limits::lowest(); - double mincoeff = (std::numeric_limits::max)(); - for (auto i = 0; i < R.size(); ++i) - { - if (std::isnan(R(i))) - continue; - if (R(i) > maxcoeff) - maxcoeff = R(i); - if (R(i) < mincoeff) - mincoeff = R(i); - sum += R(i); - } - - log()->get(LogLevel::Debug) << "R: max=" << maxcoeff - << "; min=" << mincoeff - << "; sum=" << sum - << "; size=" << R.size() << std::endl; - - // median takes an unsorted vector, possibly containing NANs, and - // returns the median value. - auto median = [&](std::vector vals) - { - // Begin by partitioning the vector by isnan. - auto ptr = std::partition(vals.begin(), vals.end(), [](double p) - { - return std::isnan(p); - }); - - // Copy the actual values, thus eliminating NANs, and sort it. - std::vector cp(ptr, vals.end()); - std::sort(cp.begin(), cp.end()); - - std::cerr << "median troubleshooting\n"; - std::cerr << vals.size() << "\t" << cp.size() << std::endl; - std::cerr << cp.size() % 2 << std::endl; - std::cerr << cp[cp.size()/2-1] << "\t" << - cp[cp.size()/2] << std::endl; - if (l == 7) - { - for (auto const& v : cp) - std::cerr << v << ", "; - std::cerr << std::endl; - } - - // Compute the median value. For even sized vectors, this is the - // average of the midpoints, otherwise it is the midpoint. - double median = 0.0; - if (cp.size() % 2 == 0) - median = (cp[cp.size()/2-1]+cp[cp.size()/2])/2; - else - median = cp[cp.size()/2]; - - return median; - }; - - // Compute median of residuals. - std::vector allres(R.data(), R.data()+R.size()); - double m = median(allres); - - // Compute absolute difference of the residuals from the median. - ArrayXd ad = (R.array()-m).abs(); - - // Compute median of absolute differences, with scale factor (1.4862) - // for a normal distribution. - std::vector absdiff(ad.data(), ad.data()+ad.size()); - double mad = 1.4862 * median(absdiff); - - // Divide absolute differences by MAD. Values greater than 2 are - // considered outliers. - MatrixXd M = (ad / mad).matrix(); - - sum = 0.0; - maxcoeff = std::numeric_limits::lowest(); - mincoeff = (std::numeric_limits::max)(); - for (auto i = 0; i < M.size(); ++i) - { - if (std::isnan(M(i))) - continue; - if (M(i) > maxcoeff) - maxcoeff = M(i); - if (M(i) < mincoeff) - mincoeff = M(i); - sum += M(i); - } - - log()->get(LogLevel::Debug) << "M: max=" << maxcoeff - << "; min=" << mincoeff - << "; sum=" << sum - << "; size=" << M.size() << std::endl; - - double madthresh = 2.0; - // Just computing the percent outlier FYI. - double perc = static_cast((M.array() > madthresh).count()); - perc /= static_cast(R.size()); - perc *= 100.0; - log()->get(LogLevel::Debug) << "median=" << m - << "; MAD=" << mad - << "; " << (M.array() > madthresh).count() - << " outliers out of " << R.size() - << " control points (" << perc << "%)\n"; - - // If the TPS control-point is recognized as a non-ground point, it is - // replaced by the interpolated point. The time complexity of the - // approach is reduced by filtering only the control-points in each - // iteration. - if (l < 3) - { - for (auto i = 0; i < M.size(); ++i) - { - if (M(i) > madthresh) - z_samp(i) = std::numeric_limits::quiet_NaN(); - // z_samp(i) = surface(i); - } - } - - if (log()->getLevel() > LogLevel::Debug5) - { - char buffer[256]; - sprintf(buffer, "interp_surface_%d.laz", l); - std::string name(buffer); - // eigen::writeMatrix(surface, name, cur_cell_size, m_bounds, srs); - writeControl(x_samp, y_samp, surface, name); - - char bufm[256]; - sprintf(bufm, "master_control_%d.laz", l); - std::string namem(bufm); - writeControl(cx, cy, cz, namem); - - // this is identical to filtered control when written here - should move it... - char buf3[256]; - sprintf(buf3, "prev_control_%d.laz", l); - std::string name3(buf3); - writeControl(x_prev, y_prev, z_prev, name3); - - char rbuf[256]; - sprintf(rbuf, "residual_%d.laz", l); - std::string rbufn(rbuf); - // eigen::writeMatrix(R, rbufn, cur_cell_size, m_bounds, srs); - writeControl(x_samp, y_samp, R, rbufn); - - char mbuf[256]; - sprintf(mbuf, "median_%d.laz", l); - std::string mbufn(mbuf); - // eigen::writeMatrix(M, mbufn, cur_cell_size, m_bounds, srs); - writeControl(x_samp, y_samp, M, mbufn); - - char buf2[256]; - sprintf(buf2, "adjusted_control_%d.laz", l); - std::string name2(buf2); - writeControl(x_samp, y_samp, z_samp, name2); - } - - x_prev = x_samp; - y_prev = y_samp; - z_prev = z_samp; - } - - MatrixXd surface = eigen::computeSpline(x_prev, y_prev, z_prev, cx, cy); - - if (log()->getLevel() > LogLevel::Debug5) - { - // writeControl(cx, cy, mc, "closed.laz"); - // - char buffer[256]; - sprintf(buffer, "final_surface.tif"); - std::string name(buffer); - eigen::writeMatrix(surface, name, "GTiff", m_cellSize, m_bounds, srs); - // - // char rbuf[256]; - // sprintf(rbuf, "final_residual.tif"); - // std::string rbufn(rbuf); - // eigen::writeMatrix(R, rbufn, cur_cell_size, m_bounds, srs); - // - // char obuf[256]; - // sprintf(obuf, "final_opened.tif"); - // std::string obufn(obuf); - // eigen::writeMatrix(maxZ, obufn, cur_cell_size, m_bounds, srs); - // - // char Tbuf[256]; - // sprintf(Tbuf, "final_tophat.tif"); - // std::string Tbufn(Tbuf); - // eigen::writeMatrix(T, Tbufn, cur_cell_size, m_bounds, srs); - // - // char tbuf[256]; - // sprintf(tbuf, "final_thresh.tif"); - // std::string tbufn(tbuf); - // eigen::writeMatrix(t, tbufn, cur_cell_size, m_bounds, srs); - } - - // apply final filtering (top hat) using raw points against TPS - - // ...the LiDAR points are filtered only at the bottom level. - for (point_count_t i = 0; i < np; ++i) - { - using namespace Dimension; - - double x = view->getFieldAs(Id::X, i); - double y = view->getFieldAs(Id::Y, i); - double z = view->getFieldAs(Id::Z, i); - - int c = Utils::clamp(getColIndex(x, cur_cell_size), 0, m_numCols-1); - int r = Utils::clamp(getRowIndex(y, cur_cell_size), 0, m_numRows-1); - - double res = z - surface(r, c); - if (res < 1.0) - groundIdx.push_back(i); - } - - return groundIdx; -} - -void MongusFilter::downsampleMin(Eigen::MatrixXd *cx, Eigen::MatrixXd *cy, - Eigen::MatrixXd* cz, Eigen::MatrixXd *dcx, - Eigen::MatrixXd *dcy, Eigen::MatrixXd* dcz, - double cell_size) -{ - int nr = static_cast(std::ceil(cz->rows() / cell_size)); - int nc = static_cast(std::ceil(cz->cols() / cell_size)); - - // std::cerr << nr << "\t" << nc << "\t" << cell_size << std::endl; - - dcx->resize(nr, nc); - dcx->setConstant(std::numeric_limits::quiet_NaN()); - - dcy->resize(nr, nc); - dcy->setConstant(std::numeric_limits::quiet_NaN()); - - dcz->resize(nr, nc); - dcz->setConstant((std::numeric_limits::max)()); - - for (auto c = 0; c < cz->cols(); ++c) - { - for (auto r = 0; r < cz->rows(); ++r) - { - if ((*cz)(r, c) == (std::numeric_limits::max)()) - continue; - - int rr = static_cast(std::floor(r / cell_size)); - int cc = static_cast(std::floor(c / cell_size)); - - if ((*cz)(r, c) < (*dcz)(rr, cc)) - { - (*dcx)(rr, cc) = (*cx)(r, c); - (*dcy)(rr, cc) = (*cy)(r, c); - (*dcz)(rr, cc) = (*cz)(r, c); - } - } - } -} - -PointViewSet MongusFilter::run(PointViewPtr view) -{ - bool logOutput = log()->getLevel() > LogLevel::Debug1; - if (logOutput) - log()->floatPrecision(8); - log()->get(LogLevel::Debug2) << "Process MongusFilter...\n"; - - std::vector idx = processGround(view); - std::cerr << idx.size() << std::endl; - - PointViewSet viewSet; - - if (!idx.empty() && (m_classify || m_extract)) - { - - if (m_classify) - { - log()->get(LogLevel::Debug2) << "Labeled " << idx.size() << - " ground returns!\n"; - - // set the classification label of ground returns as 2 - // (corresponding to ASPRS LAS specification) - for (const auto& i : idx) - { - view->setField(Dimension::Id::Classification, i, 2); - } - - viewSet.insert(view); - } - - if (m_extract) - { - log()->get(LogLevel::Debug2) << "Extracted " << idx.size() << - " ground returns!\n"; - - // create new PointView containing only ground returns - PointViewPtr output = view->makeNew(); - for (const auto& i : idx) - { - output->appendPoint(*view, i); - } - - viewSet.erase(view); - viewSet.insert(output); - } - } - else - { - if (idx.empty()) - log()->get(LogLevel::Debug2) << "Filtered cloud has no " - "ground returns!\n"; - - if (!(m_classify || m_extract)) - log()->get(LogLevel::Debug2) << "Must choose --classify or " - "--extract\n"; - - // return the view buffer unchanged - viewSet.insert(view); - } - - return viewSet; -} - -} // namespace pdal diff --git a/filters/MongusFilter.hpp b/filters/MongusFilter.hpp deleted file mode 100644 index 9473eceebf..0000000000 --- a/filters/MongusFilter.hpp +++ /dev/null @@ -1,89 +0,0 @@ -/****************************************************************************** -* Copyright (c) 2016, Bradley J Chambers (brad.chambers@gmail.com) -* -* All rights reserved. -* -* Redistribution and use in source and binary forms, with or without -* modification, are permitted provided that the following -* conditions are met: -* -* * Redistributions of source code must retain the above copyright -* notice, this list of conditions and the following disclaimer. -* * Redistributions in binary form must reproduce the above copyright -* notice, this list of conditions and the following disclaimer in -* the documentation and/or other materials provided -* with the distribution. -* * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the -* names of its contributors may be used to endorse or promote -* products derived from this software without specific prior -* written permission. -* -* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, -* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, -* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS -* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED -* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT -* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY -* OF SUCH DAMAGE. -****************************************************************************/ - -#pragma once - -#include - -#include - -#include -#include - -namespace pdal -{ - -class PointLayout; -class PointView; - -typedef std::unordered_map> PointIdHash; - -// NOTE: This filter is deprecated as of 3/2019. - -class PDAL_DLL MongusFilter : public Filter -{ -public: - MongusFilter() : Filter() - {} - - std::string getName() const; - -private: - bool m_classify; - bool m_extract; - int m_numRows; - int m_numCols; - int m_maxRow; - double m_cellSize; - double m_k; - int m_l; - BOX2D m_bounds; - - virtual void addDimensions(PointLayoutPtr layout); - virtual void addArgs(ProgramArgs& args); - int getColIndex(double x, double cell_size); - int getRowIndex(double y, double cell_size); - void writeControl(Eigen::MatrixXd cx, Eigen::MatrixXd cy, Eigen::MatrixXd cz, std::string filename); - void downsampleMin(Eigen::MatrixXd *cx, Eigen::MatrixXd *cy, - Eigen::MatrixXd* cz, Eigen::MatrixXd *dcx, - Eigen::MatrixXd *dcy, Eigen::MatrixXd* dcz, - double cell_size); - std::vector processGround(PointViewPtr view); - virtual PointViewSet run(PointViewPtr view); - - MongusFilter& operator=(const MongusFilter&); // not implemented - MongusFilter(const MongusFilter&); // not implemented -}; - -} // namespace pdal From 7f49aaefed883a380b1a0cee8bee8235fd61c94b Mon Sep 17 00:00:00 2001 From: chambbj Date: Thu, 16 May 2019 09:45:15 -0500 Subject: [PATCH 3/6] Remove KDistance filter and references, replaced by NNDistance (#2518) --- doc/stages/filters.kdistance.rst | 40 -------------- doc/stages/filters.rst | 5 -- filters/KDistanceFilter.cpp | 92 -------------------------------- filters/KDistanceFilter.hpp | 68 ----------------------- 4 files changed, 205 deletions(-) delete mode 100644 doc/stages/filters.kdistance.rst delete mode 100644 filters/KDistanceFilter.cpp delete mode 100644 filters/KDistanceFilter.hpp diff --git a/doc/stages/filters.kdistance.rst b/doc/stages/filters.kdistance.rst deleted file mode 100644 index 98b6469f8a..0000000000 --- a/doc/stages/filters.kdistance.rst +++ /dev/null @@ -1,40 +0,0 @@ -.. _filters.kdistance: - -=============================================================================== -filters.kdistance -=============================================================================== - -The **K-Distance filter** creates a new attribute ``KDistance`` that -contains the euclidean distance to a point's k-th nearest neighbor. - -.. note:: - - The K-distance filter is deprecated and has been replaced by - :ref:`filters.nndistance`. - -.. embed:: - -Example -------------------------------------------------------------------------------- - -.. code-block:: json - - [ - "input.las", - { - "type":"filters.kdistance", - "k":8 - }, - { - "type":"writers.bpf", - "filename":"output.las", - "output_dims":"X,Y,Z,KDistance" - } - ] - -Options -------------------------------------------------------------------------------- - -k - The number of k nearest neighbors. [Default: 10] - diff --git a/doc/stages/filters.rst b/doc/stages/filters.rst index 4133f10102..812f4d4b85 100644 --- a/doc/stages/filters.rst +++ b/doc/stages/filters.rst @@ -43,7 +43,6 @@ invalidate an existing KD-tree. filters.ferry filters.hag filters.info - filters.kdistance filters.lof filters.neighborclassifier filters.nndistance @@ -89,10 +88,6 @@ invalidate an existing KD-tree. Compute pointwise height above ground estimate. Requires points to be classified as ground/non-ground prior to estimating. -:ref:`filters.kdistance` - Compute pointwise K-Distance (the Euclidean distance to a point's k-th - nearest neighbor). [Deprecated - use :ref:`filters.nndistance`] - :ref:`filters.lof` Compute pointwise Local Outlier Factor (along with K-Distance and Local Reachability Distance). diff --git a/filters/KDistanceFilter.cpp b/filters/KDistanceFilter.cpp deleted file mode 100644 index 3dbd052b38..0000000000 --- a/filters/KDistanceFilter.cpp +++ /dev/null @@ -1,92 +0,0 @@ -/****************************************************************************** -* Copyright (c) 2016, Bradley J Chambers (brad.chambers@gmail.com) -* -* All rights reserved. -* -* Redistribution and use in source and binary forms, with or without -* modification, are permitted provided that the following -* conditions are met: -* -* * Redistributions of source code must retain the above copyright -* notice, this list of conditions and the following disclaimer. -* * Redistributions in binary form must reproduce the above copyright -* notice, this list of conditions and the following disclaimer in -* the documentation and/or other materials provided -* with the distribution. -* * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the -* names of its contributors may be used to endorse or promote -* products derived from this software without specific prior -* written permission. -* -* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, -* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, -* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS -* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED -* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT -* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY -* OF SUCH DAMAGE. -****************************************************************************/ - -#include "KDistanceFilter.hpp" - -#include - -#include -#include - -namespace pdal -{ - -static PluginInfo const s_info -{ - "filters.kdistance", - "K-Distance Filter", - "http://pdal.io/stages/filters.kdistance.html" -}; - -CREATE_STATIC_STAGE(KDistanceFilter, s_info) - -std::string KDistanceFilter::getName() const -{ - return s_info.name; -} - -void KDistanceFilter::addArgs(ProgramArgs& args) -{ - args.add("k", "k neighbors", m_k, 10); -} - -void KDistanceFilter::addDimensions(PointLayoutPtr layout) -{ - using namespace Dimension; - m_kdist = layout->registerOrAssignDim("KDistance", Type::Double); -} - -void KDistanceFilter::filter(PointView& view) -{ - using namespace Dimension; - - KD3Index& index = view.build3dIndex(); - - // Increment the minimum number of points, as knnSearch will be returning - // the neighbors along with the query point. - m_k++; - - // Compute the k-distance for each point. The k-distance is the Euclidean - // distance to k-th nearest neighbor. - log()->get(LogLevel::Debug) << "Computing k-distances...\n"; - for (PointId i = 0; i < view.size(); ++i) - { - std::vector indices(m_k); - std::vector sqr_dists(m_k); - index.knnSearch(i, m_k, &indices, &sqr_dists); - view.setField(m_kdist, i, std::sqrt(sqr_dists[m_k-1])); - } -} - -} // namespace pdal diff --git a/filters/KDistanceFilter.hpp b/filters/KDistanceFilter.hpp deleted file mode 100644 index 37e49e3419..0000000000 --- a/filters/KDistanceFilter.hpp +++ /dev/null @@ -1,68 +0,0 @@ -/****************************************************************************** -* Copyright (c) 2016, Bradley J Chambers (brad.chambers@gmail.com) -* -* All rights reserved. -* -* Redistribution and use in source and binary forms, with or without -* modification, are permitted provided that the following -* conditions are met: -* -* * Redistributions of source code must retain the above copyright -* notice, this list of conditions and the following disclaimer. -* * Redistributions in binary form must reproduce the above copyright -* notice, this list of conditions and the following disclaimer in -* the documentation and/or other materials provided -* with the distribution. -* * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the -* names of its contributors may be used to endorse or promote -* products derived from this software without specific prior -* written permission. -* -* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE -* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, -* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, -* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS -* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED -* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT -* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY -* OF SUCH DAMAGE. -****************************************************************************/ - -#pragma once - -#include - -#include - -namespace pdal -{ - -class PointLayout; -class PointView; -class ProgramArgs; - -class PDAL_DLL KDistanceFilter : public Filter -{ -public: - KDistanceFilter() : Filter() - {} - - std::string getName() const; - -private: - Dimension::Id m_kdist; - int m_k; - - virtual void addArgs(ProgramArgs& args); - virtual void addDimensions(PointLayoutPtr layout); - virtual void filter(PointView& view); - - KDistanceFilter& operator=(const KDistanceFilter&); // not implemented - KDistanceFilter(const KDistanceFilter&); // not implemented -}; - -} // namespace pdal From 416f99516221925d250956f30e57774b90ef5539 Mon Sep 17 00:00:00 2001 From: chambbj Date: Thu, 16 May 2019 09:46:55 -0500 Subject: [PATCH 4/6] Shift to using double precision Eigen vectors and matrices (#2519) --- filters/ApproximateCoplanarFilter.cpp | 2 +- filters/CovarianceFeaturesFilter.cpp | 28 +++++++++++++-------------- filters/EigenvaluesFilter.cpp | 2 +- filters/GreedyProjection.cpp | 2 +- filters/NormalFilter.cpp | 6 +++--- pdal/EigenUtils.cpp | 8 ++++---- pdal/EigenUtils.hpp | 2 +- 7 files changed, 25 insertions(+), 25 deletions(-) diff --git a/filters/ApproximateCoplanarFilter.cpp b/filters/ApproximateCoplanarFilter.cpp index b51961b45f..cf5d17104a 100644 --- a/filters/ApproximateCoplanarFilter.cpp +++ b/filters/ApproximateCoplanarFilter.cpp @@ -90,7 +90,7 @@ void ApproximateCoplanarFilter::filter(PointView& view) auto B = eigen::computeCovariance(view, ids); // perform the eigen decomposition - SelfAdjointEigenSolver solver(B); + SelfAdjointEigenSolver solver(B); if (solver.info() != Success) throwError("Cannot perform eigen decomposition."); auto ev = solver.eigenvalues(); diff --git a/filters/CovarianceFeaturesFilter.cpp b/filters/CovarianceFeaturesFilter.cpp index 9d79b5b517..212e35673a 100644 --- a/filters/CovarianceFeaturesFilter.cpp +++ b/filters/CovarianceFeaturesFilter.cpp @@ -77,7 +77,7 @@ void CovarianceFeaturesFilter::addDimensions(PointLayoutPtr layout) if (m_featureSet == "Dimensionality") { for (auto dim: {"Linearity", "Planarity", "Scattering", "Verticality"}) - m_extraDims[dim] = layout->registerOrAssignDim(dim, Dimension::Type::Float); + m_extraDims[dim] = layout->registerOrAssignDim(dim, Dimension::Type::Double); } } @@ -113,21 +113,21 @@ void CovarianceFeaturesFilter::setDimensionality(PointView &view, const PointId auto B = eigen::computeCovariance(view, ids); // perform the eigen decomposition - SelfAdjointEigenSolver solver(B); + SelfAdjointEigenSolver solver(B); if (solver.info() != Success) throwError("Cannot perform eigen decomposition."); // Extract eigenvalues and eigenvectors in decreasing order (largest eigenvalue first) auto ev = solver.eigenvalues(); - std::vector lambda = {(std::max(ev[2],0.f)), - (std::max(ev[1],0.f)), - (std::max(ev[0],0.f))}; + std::vector lambda = {(std::max(ev[2],0.0)), + (std::max(ev[1],0.0)), + (std::max(ev[0],0.0))}; if (lambda[0] == 0) throwError("Eigenvalues are all 0. Can't compute local features."); auto eigenVectors = solver.eigenvectors(); - std::vector v1(3), v2(3), v3(3); + std::vector v1(3), v2(3), v3(3); for (int i=0; i < 3; i++) { v1[i] = eigenVectors.col(2)(i); @@ -135,21 +135,21 @@ void CovarianceFeaturesFilter::setDimensionality(PointView &view, const PointId v3[i] = eigenVectors.col(0)(i); } - float linearity = (sqrtf(lambda[0]) - sqrtf(lambda[1])) / sqrtf(lambda[0]); - float planarity = (sqrtf(lambda[1]) - sqrtf(lambda[2])) / sqrtf(lambda[0]); - float scattering = sqrtf(lambda[2]) / sqrtf(lambda[0]); + double linearity = (sqrt(lambda[0]) - sqrt(lambda[1])) / sqrt(lambda[0]); + double planarity = (sqrt(lambda[1]) - sqrt(lambda[2])) / sqrt(lambda[0]); + double scattering = sqrt(lambda[2]) / sqrt(lambda[0]); view.setField(m_extraDims["Linearity"], id, linearity); view.setField(m_extraDims["Planarity"], id, planarity); view.setField(m_extraDims["Scattering"], id, scattering); - std::vector unary_vector(3); - float norm = 0; + std::vector unary_vector(3); + double norm = 0; for (int i=0; i <3 ; i++) { - unary_vector[i] = lambda[0] * fabsf(v1[i]) + lambda[1] * fabsf(v2[i]) + lambda[2] * fabsf(v3[i]); + unary_vector[i] = lambda[0] * fabs(v1[i]) + lambda[1] * fabs(v2[i]) + lambda[2] * fabs(v3[i]); norm += unary_vector[i] * unary_vector[i]; } - norm = sqrtf(norm); + norm = sqrt(norm); view.setField(m_extraDims["Verticality"], id, unary_vector[2] / norm); } -} \ No newline at end of file +} diff --git a/filters/EigenvaluesFilter.cpp b/filters/EigenvaluesFilter.cpp index edd1c0a009..535a579de0 100644 --- a/filters/EigenvaluesFilter.cpp +++ b/filters/EigenvaluesFilter.cpp @@ -89,7 +89,7 @@ void EigenvaluesFilter::filter(PointView& view) auto B = eigen::computeCovariance(view, ids); // perform the eigen decomposition - SelfAdjointEigenSolver solver(B); + SelfAdjointEigenSolver solver(B); if (solver.info() != Success) throwError("Cannot perform eigen decomposition."); auto ev = solver.eigenvalues(); diff --git a/filters/GreedyProjection.cpp b/filters/GreedyProjection.cpp index 56b3ed50d7..c6121172ec 100644 --- a/filters/GreedyProjection.cpp +++ b/filters/GreedyProjection.cpp @@ -188,7 +188,7 @@ void GreedyProjection::filter(PointView& view) std::vector point2index (input_->points.size (), -1); for (int cp = 0; cp < static_cast (indices_->size ()); ++cp) { - coords_.push_back(input_->points[(*indices_)[cp]].getVector3fMap()); + coords_.push_back(input_->points[(*indices_)[cp]].getVector3dMap()); point2index[(*indices_)[cp]] = cp; } **/ diff --git a/filters/NormalFilter.cpp b/filters/NormalFilter.cpp index ce4bfe946f..7b4e7c52d7 100644 --- a/filters/NormalFilter.cpp +++ b/filters/NormalFilter.cpp @@ -127,18 +127,18 @@ void NormalFilter::filter(PointView& view) auto B = eigen::computeCovariance(view, ids); // perform the eigen decomposition - Eigen::SelfAdjointEigenSolver solver(B); + Eigen::SelfAdjointEigenSolver solver(B); if (solver.info() != Eigen::Success) throwError("Cannot perform eigen decomposition."); auto eval = solver.eigenvalues(); - Eigen::Vector3f normal = solver.eigenvectors().col(0); + Eigen::Vector3d normal = solver.eigenvectors().col(0); if (m_viewpointArg->set()) { using namespace Dimension; PointRef p = view.point(i); - Eigen::Vector3f vp( + Eigen::Vector3d vp( (float)(m_args->m_viewpoint.x() - p.getFieldAs(Id::X)), (float)(m_args->m_viewpoint.y() - p.getFieldAs(Id::Y)), (float)(m_args->m_viewpoint.z() - p.getFieldAs(Id::Z))); diff --git a/pdal/EigenUtils.cpp b/pdal/EigenUtils.cpp index 989ec7e511..a8f4a294d0 100644 --- a/pdal/EigenUtils.cpp +++ b/pdal/EigenUtils.cpp @@ -83,7 +83,7 @@ Eigen::Vector3d computeCentroid(PointView& view, return centroid; } -Eigen::Matrix3f computeCovariance(PointView& view, +Eigen::Matrix3d computeCovariance(PointView& view, const std::vector& ids) { using namespace Eigen; @@ -93,7 +93,7 @@ Eigen::Matrix3f computeCovariance(PointView& view, Vector3d centroid = computeCentroid(view, ids); // demean the neighborhood - MatrixXf A(3, n); + MatrixXd A(3, n); size_t k = 0; for (auto const& j : ids) { @@ -117,9 +117,9 @@ uint8_t computeRank(PointView& view, const std::vector& ids, { using namespace Eigen; - Matrix3f B = computeCovariance(view, ids); + Matrix3d B = computeCovariance(view, ids); - JacobiSVD svd(B); + JacobiSVD svd(B); svd.setThreshold((float)threshold); return static_cast(svd.rank()); diff --git a/pdal/EigenUtils.hpp b/pdal/EigenUtils.hpp index ba068f9cf3..3ca8c01a1d 100644 --- a/pdal/EigenUtils.hpp +++ b/pdal/EigenUtils.hpp @@ -100,7 +100,7 @@ PDAL_DLL Eigen::Vector3d computeCentroid(PointView& view, \param ids a vector of PointIds specifying a subset of points. \return the covariance matrix of the XYZ dimensions. */ -PDAL_DLL Eigen::Matrix3f computeCovariance(PointView& view, +PDAL_DLL Eigen::Matrix3d computeCovariance(PointView& view, const std::vector& ids); /** From 305f56c9fc5bf3b037693f450ba1d72aaf135e24 Mon Sep 17 00:00:00 2001 From: chambbj Date: Thu, 16 May 2019 09:47:33 -0500 Subject: [PATCH 5/6] Add covariancefeatures to the main filters documentation (#2517) --- doc/stages/filters.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/doc/stages/filters.rst b/doc/stages/filters.rst index 812f4d4b85..65e7ae8e6b 100644 --- a/doc/stages/filters.rst +++ b/doc/stages/filters.rst @@ -36,6 +36,7 @@ invalidate an existing KD-tree. filters.cluster filters.colorinterp filters.colorization + filters.covariancefeatures filters.dem filters.eigenvalues filters.estimaterank @@ -72,6 +73,10 @@ invalidate an existing KD-tree. :ref:`filters.colorization` Fetch and assign RGB color information from a GDAL-readable datasource. +:ref:`filters.covariancefeatures` + Filter that calculates local features based on the covariance matrix of a + point's neighborhood. + :ref:`filters.eigenvalues` Compute pointwise eigenvalues, based on k-nearest neighbors. From 31e5cca10d7fae909b5cfd93adac15bfafd4a04c Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Fri, 24 May 2019 13:01:08 -0400 Subject: [PATCH 6/6] Read data from the metadata child "data" rather than the VLR node's data for forwarded VLRs. --- io/LasReader.cpp | 2 +- io/LasWriter.cpp | 3 ++- test/unit/io/LasReaderTest.cpp | 7 ++++--- test/unit/io/LasWriterTest.cpp | 13 +++++++++++++ 4 files changed, 20 insertions(+), 5 deletions(-) diff --git a/io/LasReader.cpp b/io/LasReader.cpp index 445f41fc2a..97198cdbab 100644 --- a/io/LasReader.cpp +++ b/io/LasReader.cpp @@ -502,7 +502,6 @@ void LasReader::extractVlrMetadata(MetadataNode& forward, MetadataNode& m) std::ostringstream name; name << "vlr_" << i++; MetadataNode vlrNode(name.str()); - m.add(vlrNode); vlrNode.addEncoded("data", (const uint8_t *)vlr.data(), vlr.dataLen(), vlr.description()); @@ -512,6 +511,7 @@ void LasReader::extractVlrMetadata(MetadataNode& forward, MetadataNode& m) vlrNode.add("record_id", vlr.recordId(), "Record ID specified by the user."); vlrNode.add("description", vlr.description()); + m.add(vlrNode); if (vlr.userId() == TRANSFORM_USER_ID|| vlr.userId() == LASZIP_USER_ID || diff --git a/io/LasWriter.cpp b/io/LasWriter.cpp index 763f9b1f72..5880a09f83 100644 --- a/io/LasWriter.cpp +++ b/io/LasWriter.cpp @@ -466,7 +466,8 @@ void LasWriter::addForwardVlrs() const MetadataNode& recordIdNode = n.findChild("record_id"); if (recordIdNode.valid() && userIdNode.valid()) { - data = Utils::base64_decode(n.value()); + const MetadataNode& dataNode = n.findChild("data"); + data = Utils::base64_decode(dataNode.value()); uint16_t recordId = (uint16_t)std::stoi(recordIdNode.value()); addVlr(userIdNode.value(), recordId, n.description(), data); } diff --git a/test/unit/io/LasReaderTest.cpp b/test/unit/io/LasReaderTest.cpp index 7bcba1bf3b..af57215a33 100644 --- a/test/unit/io/LasReaderTest.cpp +++ b/test/unit/io/LasReaderTest.cpp @@ -198,9 +198,11 @@ TEST(LasReaderTest, inspect) QuickInfo qi = reader.preview(); - std::string testWkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]"; + std::string testWkt { + R"(GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]])" + }; - EXPECT_EQ(qi.m_srs.getWKT(), testWkt); + EXPECT_TRUE(Utils::startsWith(qi.m_srs.getWKT(), testWkt)); EXPECT_EQ(qi.m_pointCount, 5380u); BOX3D bounds(-94.683465399999989, 31.0367341, 39.081000199999998, @@ -531,5 +533,4 @@ TEST(LasReaderTest, IgnoreVLRs) m = m.findChild("data"); EXPECT_FALSE(!m.empty()) << "No value for node " << i; } - } diff --git a/test/unit/io/LasWriterTest.cpp b/test/unit/io/LasWriterTest.cpp index 727aede242..21e7ea97a0 100644 --- a/test/unit/io/LasWriterTest.cpp +++ b/test/unit/io/LasWriterTest.cpp @@ -1235,9 +1235,22 @@ TEST(LasWriterTest, forward_spec_3) !temp.findChild(recPred).empty() && !temp.findChild(userPred).empty(); }; + MetadataNode origRoot = reader.getMetadata(); + MetadataNodeList origNodes = origRoot.findChildren(pred); + EXPECT_EQ(origNodes.size(), 1u); + MetadataNode origNode = origNodes[0]; + MetadataNode root = reader2.getMetadata(); MetadataNodeList nodes = root.findChildren(pred); EXPECT_EQ(nodes.size(), 1u); + MetadataNode node = nodes[0]; + + // Also test that we're properly forwarding data. + origNode = origNode.findChild("data"); + node = node.findChild("data"); + EXPECT_EQ(origNode.value().size(), 28u); + EXPECT_EQ(node.value().size(), origNode.value().size()); + EXPECT_EQ(node.value(), origNode.value()); } TEST(LasWriterTest, oversize_vlr)