diff --git a/doc/references.rst b/doc/references.rst index c894fa1dd0..43c2249f23 100644 --- a/doc/references.rst +++ b/doc/references.rst @@ -38,8 +38,12 @@ Reference .. [Cook1986] Cook, Robert L. "Stochastic sampling in computer graphics." *ACM Transactions on Graphics (TOG)* 5.1 (1986): 51-72. +.. [Demantke2011] Demantké J., Mallet C., David N., Vallet, B. "Dimensionality Based Scale Selection in 3d LIDAR Point Clouds." Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci, XXXVIII-5/W12, 97-102, 2011 + .. [Dippe1985] Dippé, Mark AZ, and Erling Henry Wold. "Antialiasing through stochastic sampling." *ACM Siggraph Computer Graphics* 19.3 (1985): 69-78. +.. [Guinard2017] Guinard S., Landrieu L. "Weakly Supervised Segmented-Aided Classification of Urban Scenes From 3D LIDAR Point Clouds." Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., XLII-1/W1, 151-157, 2017 + .. [Kazhdan2006] Kazhdan, Michael, Matthew Bolitho, and Hugues Hoppe. "Poisson surface reconstruction." Proceedings of the fourth Eurographics symposium on Geometry processing. Vol. 7. 2006. .. [Li2010] Li, Ruosi, et al. "Polygonizing extremal surfaces with manifold guarantees." Proceedings of the 14th ACM Symposium on Solid and Physical Modeling. ACM, 2010. diff --git a/doc/stages/filters.covariancefeatures.rst b/doc/stages/filters.covariancefeatures.rst new file mode 100644 index 0000000000..1c1d215c31 --- /dev/null +++ b/doc/stages/filters.covariancefeatures.rst @@ -0,0 +1,65 @@ +.. _filters.covariancefeatures: + +=============================================================================== +filters.covariancefeatures +=============================================================================== + +This filter implements various local feature descriptors introduced that are based on the covariance matrix of a +point's neighborhood. The user can pick a set of feature descriptors by +setting the ``feature_set`` option. Currently, the only supported feature is the dimensionality_ +set of feature descriptors introduced below. + +Example +------------------------------------------------------------------------------- + +.. code-block:: json + + [ + "input.las", + { + "type":"filters.covariancefeatures", + "knn":8, + "threads": 2, + "feature_set": "Dimensionality" + }, + { + "type":"writers.bpf", + "filename":"output.las", + "output_dims":"X,Y,Z,Linearity,Planarity,Scattering,Verticality" + } + ] + +Options +------------------------------------------------------------------------------- + +knn + The number of k nearest neighbors used for calculating the covariance matrix. [Default: 10] + +threads + The number of threads used for computing the feature descriptors. [Default: 1] + +feature_set + The features to be computed. Currently only supports ``Dimensionality``. [Default: "Dimensionality"] + +.. _dimensionality: + +Dimensionality feature set +................................................................................ +The features introduced in [Demantke2011]_ describe the shape +of the neighborhood, indicating whether +the local geometry is more linear (1D), planar (2D) or volumetric (3D) while the one introduced in +[Guinard2017]_ adds the idea of a structure being vertical. + + +The dimensionality filter introduces the following four descriptors that are computed from the covariance matrix of the ``knn`` neighbors: + +* linearity - higher for long thin strips +* planarity - higher for planar surfaces +* scattering - higher for complex 3d neighbourhoods +* verticality - higher for vertical structures, highest for thin vertical strips + +It introduces four new dimensions that hold each one of these values: ``Linearity`` ``Planarity`` ``Scattering`` +and ``Verticality``. + + + diff --git a/doc/stages/filters.hexbin.rst b/doc/stages/filters.hexbin.rst index f6e8d73c41..0f298c9806 100644 --- a/doc/stages/filters.hexbin.rst +++ b/doc/stages/filters.hexbin.rst @@ -115,5 +115,5 @@ threshold is considered "in" the data set. [Default: 15] precision - Coordinate precision to use in writing out the well-known text - of the boundary polygon. [Default: 8] + Minimum number of significant digits to use in writing out the + well-known text of the boundary polygon. [Default: 8] diff --git a/doc/stages/filters.reprojection.rst b/doc/stages/filters.reprojection.rst index 59d03fe381..3bdd92e0c2 100644 --- a/doc/stages/filters.reprojection.rst +++ b/doc/stages/filters.reprojection.rst @@ -24,7 +24,7 @@ reprojecting. Example 1 -------------------------------------------------------------------------------- -This pipeline reprojecs terrain points with Z-values between 0 and 100 by first +This pipeline reprojects terrain points with Z-values between 0 and 100 by first applying a range filter and then specifing both the input and output spatial reference as EPSG-codes. The X and Y dimensions are scaled to allow enough precision in the output coordinates. diff --git a/doc/stages/writers.ply.rst b/doc/stages/writers.ply.rst index 8bb957087e..012e8a7bf3 100644 --- a/doc/stages/writers.ply.rst +++ b/doc/stages/writers.ply.rst @@ -42,12 +42,21 @@ storage_mode of the machine. [Default: "ascii"] dims - List of dimensions to write as elements. [Default: all dimensions] + List of dimensions (and sizes) in the format [=],... + to write as output. (e.g., "Y=int32_t, X,Red=char") + [Default: All dimensions with stored types] faces Write a mesh as faces in addition to writing points as vertices. [Default: false] +sized_types + PLY has variously been written with explicitly sized type strings + ('int8', 'float32", 'uint32', etc.) and implied sized type strings + ('char', 'float', 'int', etc.). If true, explicitly sized type strings + are used. If false, implicitly sized type strings are used. + [Default: true] + 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'. diff --git a/doc/tutorial/iowa-entwine.rst b/doc/tutorial/iowa-entwine.rst index 75ae8e052b..e08d69e4c1 100644 --- a/doc/tutorial/iowa-entwine.rst +++ b/doc/tutorial/iowa-entwine.rst @@ -10,7 +10,7 @@ This tutorial describes how to use `Conda`_, `Entwine`_, `PDAL`_, and `GDAL`_ to read data from the `USGS 3DEP AWS Public Dataset`_. We will be using PDAL's `readers.ept`_ to fetch data, we will filter it for noise using `filters.outlier`_, we will classify the data as ground/not-ground using `filters.smrf`_, and we will -write out a digital terrain model with `writers.gdal`. Once our elevation model +write out a digital terrain model with :ref:`writers.gdal`. Once our elevation model is constructed, we will use GDAL `gdaldem`_ operations to create hillshade, slope, and color relief. @@ -150,7 +150,7 @@ each point to ``0``. .. figure:: ../images/pipeline-example-filters.assign.png :scale: 50% - :ref`filters.assign` can also take in an option to apply assignments + :ref:`filters.assign` can also take in an option to apply assignments based on a conditional. If you want to assign values based on a bounding geometry, use :ref:`filters.overlay`. diff --git a/filters/CovarianceFeaturesFilter.cpp b/filters/CovarianceFeaturesFilter.cpp new file mode 100644 index 0000000000..9d79b5b517 --- /dev/null +++ b/filters/CovarianceFeaturesFilter.cpp @@ -0,0 +1,155 @@ +/****************************************************************************** +* Copyright (c) 2019, Helix Re Inc. nicolas@helix.re +* +* 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 Helix Re Inc. 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. +****************************************************************************/ + +// This is an implementation of the local feature descriptors introduced in +// WEAKLY SUPERVISED SEGMENTATION-AIDED CLASSIFICATION OF URBANSCENES FROM 3D LIDAR POINT CLOUDS +// Stéphane Guinard, Loïc Landrieu, 2017 + +#include "CovarianceFeaturesFilter.hpp" + +#include +#include +#include + +#include + +#include +#include +#include + +namespace pdal +{ + +static StaticPluginInfo const s_info +{ + "filters.covariancefeatures", + "Filter that calculates local features based on the covariance matrix of a point's neighborhood.", + "http://pdal.io/stages/filters.covariancefeatures.html" +}; + +CREATE_STATIC_STAGE(CovarianceFeaturesFilter, s_info) + +std::string CovarianceFeaturesFilter::getName() const +{ + return s_info.name; +} + +void CovarianceFeaturesFilter::addArgs(ProgramArgs& args) +{ + args.add("knn", "k-Nearest neighbors", m_knn, 10); + args.add("threads", "Number of threads used to run this filter", m_threads, 1); + args.add("feature_set", "Set of features to be computed", m_featureSet, "Dimensionality"); +} + +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); + } +} + +void CovarianceFeaturesFilter::filter(PointView& view) +{ + + KD3Index& kdi = view.build3dIndex(); + + point_count_t nloops = view.size(); + std::vector threadPool(m_threads); + for(int t = 0;t 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))}; + + 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); + for (int i=0; i < 3; i++) + { + v1[i] = eigenVectors.col(2)(i); + v2[i] = eigenVectors.col(1)(i); + 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]); + 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; + 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]); + norm += unary_vector[i] * unary_vector[i]; + } + norm = sqrtf(norm); + view.setField(m_extraDims["Verticality"], id, unary_vector[2] / norm); +} +} \ No newline at end of file diff --git a/filters/CovarianceFeaturesFilter.hpp b/filters/CovarianceFeaturesFilter.hpp new file mode 100644 index 0000000000..14274acaf3 --- /dev/null +++ b/filters/CovarianceFeaturesFilter.hpp @@ -0,0 +1,66 @@ +/****************************************************************************** +* Copyright (c) 2019, Helix Re Inc. nicolas@helix.re +* +* 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 Helix Re Inc. 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 PDAL_DLL CovarianceFeaturesFilter: public Filter +{ +public: + CovarianceFeaturesFilter() : Filter() {} + CovarianceFeaturesFilter &operator=(const CovarianceFeaturesFilter &) = delete; + CovarianceFeaturesFilter(const CovarianceFeaturesFilter &) = delete; + + std::string getName() const; + +private: + + int m_knn; + int m_threads; + std::string m_featureSet; + std::map m_extraDims; + + virtual void addDimensions(PointLayoutPtr layout); + virtual void addArgs(ProgramArgs &args); + virtual void filter(PointView &view); + + void setDimensionality(PointView &view, const PointId &id, const KD3Index &kid); +}; +} + diff --git a/io/PlyWriter.cpp b/io/PlyWriter.cpp index 47b4246f22..5afad2cee0 100644 --- a/io/PlyWriter.cpp +++ b/io/PlyWriter.cpp @@ -66,38 +66,58 @@ 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); + args.add("sized_types", + "Write types as size-explicit strings (e.g. 'uint32')", + m_sizedTypes, true); m_precisionArg = &args.add("precision", "Output precision", m_precision, 3); } void PlyWriter::prepared(PointTableRef table) { + PointLayoutPtr layout = table.layout(); + 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) + for (auto& dimspec : m_dimNames) { - auto id = table.layout()->findDim(name); + Dimension::Type type; + StringList parts = Utils::split2(dimspec, '='); + std::string name = parts[0]; + Utils::trim(name); + dimspec = Utils::tolower(name); + auto id = layout->findDim(name); if (id == Dimension::Id::Unknown) throwError("Unknown dimension '" + name + "' in provided " "dimension list."); - m_dims.push_back(id); + if (parts.size() == 1) + type = layout->dimType(id); + else if (parts.size() == 2) + { + Utils::trim(parts[1]); + type = Dimension::type(parts[1]); + } + else + throwError("Invalid format for dimension specification for " + "dimension '" + name + "'. " "Must be 'name[=type]'."); + m_dims.emplace_back(id, type); } } else { - m_dims = table.layout()->dims(); - for (auto dim : m_dims) - m_dimNames.push_back(Utils::tolower(table.layout()->dimName(dim))); + m_dims = layout->dimTypes(); + for (const auto& dim : m_dims) + m_dimNames.push_back(Utils::tolower(layout->dimName(dim.m_id))); } } std::string PlyWriter::getType(Dimension::Type type) const { - static std::map types = + static std::map sizedTypes = { { Dimension::Type::Signed8, "int8" }, { Dimension::Type::Unsigned8, "uint8" }, @@ -108,10 +128,21 @@ std::string PlyWriter::getType(Dimension::Type type) const { Dimension::Type::Float, "float32" }, { Dimension::Type::Double, "float64" } }; + static std::map types = + { + { Dimension::Type::Signed8, "char" }, + { Dimension::Type::Unsigned8, "uchar" }, + { Dimension::Type::Signed16, "short" }, + { Dimension::Type::Unsigned16, "ushort" }, + { Dimension::Type::Signed32, "int" }, + { Dimension::Type::Unsigned32, "uint" }, + { Dimension::Type::Float, "float" }, + { Dimension::Type::Double, "double" } + }; try { - return types.at(type); + return m_sizedTypes ? sizedTypes.at(type) : types.at(type); } catch (std::out_of_range&) { @@ -133,7 +164,7 @@ void PlyWriter::writeHeader(PointLayoutPtr layout) const for (auto dim : m_dims) { std::string name = *ni++; - std::string typeString = getType(layout->dimType(dim)); + std::string typeString = getType(dim.m_type); *m_stream << "property " << typeString << " " << name << std::endl; } if (m_faces) @@ -163,21 +194,85 @@ void PlyWriter::write(const PointViewPtr data) } +template +void writeTextVal(std::ostream& out, PointRef& point, Dimension::Id dim) +{ + T t = point.getFieldAs(dim); + out << t; +} + +// Writing int/uint can write "char" type instead of numeric values. +// Two specializations is easier than the SFINAE mess. +template<> +void writeTextVal(std::ostream& out, PointRef& point, Dimension::Id dim) +{ + int i = point.getFieldAs(dim); + out << i; +} + +template<> +void writeTextVal(std::ostream& out, PointRef& point, + Dimension::Id dim) +{ + uint32_t i = point.getFieldAs(dim); + out << i; +} + + void PlyWriter::writeValue(PointRef& point, Dimension::Id dim, Dimension::Type type) { if (m_format == Format::Ascii) { - double d = point.getFieldAs(dim); - if (m_precisionArg->set() && - Dimension::base(type) == Dimension::BaseType::Floating) + if (Dimension::base(type) == Dimension::BaseType::Floating) { - *m_stream << std::fixed; - m_stream->precision(m_precision); + if (m_precisionArg->set()) + { + *m_stream << std::fixed; + m_stream->precision(m_precision); + } + + if (type == Dimension::Type::Float) + writeTextVal(*m_stream, point, dim); + else + writeTextVal(*m_stream, point, dim); + + if (m_precisionArg->set()) + m_stream->unsetf(std::ios_base::fixed); } else - m_stream->unsetf(std::ios_base::fixed); - *m_stream << d; + { + switch (type) + { + case Dimension::Type::Unsigned8: + writeTextVal(*m_stream, point, dim); + break; + case Dimension::Type::Unsigned16: + writeTextVal(*m_stream, point, dim); + break; + case Dimension::Type::Unsigned32: + writeTextVal(*m_stream, point, dim); + break; + case Dimension::Type::Unsigned64: + writeTextVal(*m_stream, point, dim); + break; + case Dimension::Type::Signed8: + writeTextVal(*m_stream, point, dim); + break; + case Dimension::Type::Signed16: + writeTextVal(*m_stream, point, dim); + break; + case Dimension::Type::Signed32: + writeTextVal(*m_stream, point, dim); + break; + case Dimension::Type::Signed64: + writeTextVal(*m_stream, point, dim); + break; + default: + throwError("Internal error: invalid type found writing " + "output."); + } + } } else if (m_format == Format::BinaryLe) { @@ -200,8 +295,7 @@ void PlyWriter::writePoint(PointRef& point, PointLayoutPtr layout) { for (auto it = m_dims.begin(); it != m_dims.end();) { - Dimension::Id dim = *it; - writeValue(point, dim, layout->dimType(dim)); + writeValue(point, it->m_id, it->m_type); ++it; if (m_format == Format::Ascii && it != m_dims.end()) *m_stream << " "; diff --git a/io/PlyWriter.hpp b/io/PlyWriter.hpp index fc3e2c31b9..fc4410d51d 100644 --- a/io/PlyWriter.hpp +++ b/io/PlyWriter.hpp @@ -72,8 +72,9 @@ class PDAL_DLL PlyWriter : public Writer Format m_format; bool m_faces; StringList m_dimNames; - Dimension::IdList m_dims; + DimTypeList m_dims; int m_precision; + bool m_sizedTypes; Arg *m_precisionArg; std::vector m_views; }; diff --git a/kernels/private/density/OGR.cpp b/kernels/private/density/OGR.cpp index e65ef47fef..785d264442 100644 --- a/kernels/private/density/OGR.cpp +++ b/kernels/private/density/OGR.cpp @@ -110,7 +110,7 @@ OGRGeometryH collectHexagon(hexer::HexInfo const& info, } -} // unnamed namespace +} // unnamed namespace OGR::OGR(std::string const& filename, std::string srs, std::string driver, @@ -134,39 +134,28 @@ OGR::~OGR() void OGR::createLayer() { - OGRSFDriverH driver = OGRGetDriverByName(m_driver.c_str()); if (driver == NULL) { throw pdal::pdal_error("OGR Driver was null!"); } - char **papszDSCO = NULL; - if (pdal::FileUtils::fileExists(m_filename)) - { + if (FileUtils::fileExists(m_filename)) m_ds = OGR_Dr_Open(driver, m_filename.c_str(), TRUE /*update*/); - } else { - - m_ds = OGR_Dr_CreateDataSource(driver, m_filename.c_str(), NULL); if (m_ds == NULL) - { - throw pdal::pdal_error("Data source creation was null!"); - } + throw pdal_error("Unable to create output file '" + m_filename + + "' for density output."); } - pdal::gdal::SpatialRef spatialref(m_srs); - OGRSpatialReferenceH ref = (OGRSpatialReferenceH)spatialref.get(); - - if (!m_layerName.size()) + if (m_layerName.empty()) m_layerName = m_filename; - m_layer = OGR_DS_CreateLayer(m_ds, m_layerName.c_str(), ref, wkbMultiPolygon, NULL); + m_layer = GDALDatasetCreateLayer(m_ds, m_layerName.c_str(), m_srs.get(), + wkbMultiPolygon, NULL); if (m_layer == NULL) - { - throw pdal::pdal_error("Layer creation was null!"); - } + throw pdal_error("Layer creation was null!"); OGRFieldDefnH hFieldDefn; hFieldDefn = OGR_Fld_Create("ID", OFTInteger); @@ -188,10 +177,9 @@ void OGR::createLayer() throw pdal::pdal_error(oss.str()); } OGR_Fld_Destroy(hFieldDefn); - - } + void OGR::writeBoundary(hexer::HexGrid *grid) { OGRGeometryH multi = OGR_G_CreateGeometry(wkbMultiPolygon); diff --git a/kernels/private/density/OGR.hpp b/kernels/private/density/OGR.hpp index 12f76977f9..39ce1b2501 100644 --- a/kernels/private/density/OGR.hpp +++ b/kernels/private/density/OGR.hpp @@ -35,6 +35,8 @@ #include +#include + #include "ogr_api.h" #include "gdal.h" @@ -60,7 +62,7 @@ class OGR private: std::string m_filename; std::string m_driver; - std::string m_srs; + gdal::SpatialRef m_srs; OGRDataSourceH m_ds; OGRLayerH m_layer; diff --git a/pdal/DimUtil.hpp b/pdal/DimUtil.hpp index 9ba6788850..f639d34313 100644 --- a/pdal/DimUtil.hpp +++ b/pdal/DimUtil.hpp @@ -147,25 +147,25 @@ inline Type type(std::string s) { s = Utils::tolower(s); - if (s == "int8_t" || s == "int8") + if (s == "int8_t" || s == "int8" || s == "char") return Type::Signed8; - if (s == "int16_t" || s == "int16") + if (s == "int16_t" || s == "int16" || s == "short") return Type::Signed16; - if (s == "int32_t" || s == "int32") + if (s == "int32_t" || s == "int32" || s == "int") return Type::Signed32; - if (s == "int64_t" || s == "int64") + if (s == "int64_t" || s == "int64" || s == "long") return Type::Signed64; - if (s == "uint8_t" || s == "uint8") + if (s == "uint8_t" || s == "uint8" || s == "uchar") return Type::Unsigned8; - if (s == "uint16_t" || s == "uint16") + if (s == "uint16_t" || s == "uint16" || s == "ushort") return Type::Unsigned16; - if (s == "uint32_t" || s == "uint32") + if (s == "uint32_t" || s == "uint32" || s == "uint") return Type::Unsigned32; - if (s == "uint64_t" || s == "uint64") + if (s == "uint64_t" || s == "uint64" || s == "ulong") return Type::Unsigned64; - if (s == "float") + if (s == "float" || s == "float32") return Type::Float; - if (s == "double") + if (s == "double" || s == "float64") return Type::Double; return Type::None; } diff --git a/pdal/GDALUtils.hpp b/pdal/GDALUtils.hpp index 335ae01d39..3b8150c7c1 100644 --- a/pdal/GDALUtils.hpp +++ b/pdal/GDALUtils.hpp @@ -287,6 +287,7 @@ enum class GDALError CantOpen, NoData, InvalidBand, + BadBand, NoTransform, NotInvertible, CantReadBlock, @@ -299,6 +300,7 @@ enum class GDALError }; struct InvalidBand {}; +struct BadBand {}; struct CantReadBlock {}; struct CantWriteBlock { @@ -322,14 +324,14 @@ class Band friend class Raster; private: - GDALDataset *m_ds; /// Dataset handle - int m_bandNum; /// Band number. Band numbers start at 1. - double m_dstNoData; /// Output no data value. - GDALRasterBand *m_band; /// Band handle - int m_xTotalSize, m_yTotalSize; /// Total size (x and y) of the raster - int m_xBlockSize, m_yBlockSize; /// Size (x and y) of blocks - int m_xBlockCnt, m_yBlockCnt; /// Number of blocks in each direction - std::vector m_buf; /// Block read buffer. + GDALDataset *m_ds; /// Dataset handle + int m_bandNum; /// Band number. 1-indexed. + double m_dstNoData; /// Output no data value. + GDALRasterBand *m_band; /// Band handle + size_t m_xTotalSize, m_yTotalSize; /// Total size (x and y) of the raster + size_t m_xBlockSize, m_yBlockSize; /// Size (x and y) of blocks + size_t m_xBlockCnt, m_yBlockCnt; /// Number of blocks in each direction + std::vector m_buf; /// Block read buffer. std::string m_name; /// Band name. /** @@ -358,10 +360,19 @@ friend class Raster; m_band->SetOffset(m_band->GetOffset(NULL) - .00001); } - m_xTotalSize = m_band->GetXSize(); - m_yTotalSize = m_band->GetYSize(); + int xTotalSize = m_band->GetXSize(); + int yTotalSize = m_band->GetYSize(); - m_band->GetBlockSize(&m_xBlockSize, &m_yBlockSize); + int xBlockSize, yBlockSize; + m_band->GetBlockSize(&xBlockSize, &yBlockSize); + if (xBlockSize <= 0 || yBlockSize <= 0 || + xTotalSize <= 0 || yTotalSize <= 0) + throw BadBand(); + + m_xTotalSize = (size_t)xTotalSize; + m_yTotalSize = (size_t)yTotalSize; + m_xBlockSize = (size_t)xBlockSize; + m_yBlockSize = (size_t)yBlockSize; m_buf.resize(m_xBlockSize * m_yBlockSize); m_xBlockCnt = ((m_xTotalSize - 1) / m_xBlockSize) + 1; @@ -382,8 +393,8 @@ friend class Raster; { data.resize(m_xTotalSize * m_yTotalSize); - for (int y = 0; y < m_yBlockCnt; ++y) - for (int x = 0; x < m_xBlockCnt; ++x) + for (size_t y = 0; y < m_yBlockCnt; ++y) + for (size_t x = 0; x < m_xBlockCnt; ++x) readBlock(x, y, data); } @@ -399,19 +410,22 @@ friend class Raster; \param data Pointer to the data vector that contains the raster information. */ - void readBlock(int x, int y, std::vector& data) + void readBlock(size_t x, size_t y, std::vector& data) { uint8_t *buf = reinterpret_cast(m_buf.data()); - if (m_band->ReadBlock(x, y, buf) != CPLE_None) + + // Block indices are guaranteed not to overflow an int. + if (m_band->ReadBlock(static_cast(x), + static_cast(y), buf) != CPLE_None) throw CantReadBlock(); - int xWidth = 0; + size_t xWidth = 0; if (x == m_xBlockCnt - 1) xWidth = m_xTotalSize % m_xBlockSize; if (xWidth == 0) xWidth = m_xBlockSize; - int yHeight = 0; + size_t yHeight = 0; if (y == m_yBlockCnt - 1) yHeight = m_yTotalSize % m_yBlockSize; if (yHeight == 0) @@ -420,10 +434,10 @@ friend class Raster; auto bi = m_buf.begin(); // Go through rows copying data. Increment the buffer pointer by the // width of the row. - for (int row = 0; row < yHeight; ++row) + for (size_t row = 0; row < yHeight; ++row) { - int wholeRows = m_xTotalSize * ((y * m_yBlockSize) + row); - int partialRows = m_xBlockSize * x; + size_t wholeRows = m_xTotalSize * ((y * m_yBlockSize) + row); + size_t partialRows = m_xBlockSize * x; auto di = data.begin() + (wholeRows + partialRows); std::copy(bi, bi + xWidth, di); @@ -441,8 +455,8 @@ friend class Raster; template void write(SOURCE_ITER si, ITER_VAL srcNoData) { - for (int y = 0; y < m_yBlockCnt; ++y) - for (int x = 0; x < m_xBlockCnt; ++x) + for (size_t y = 0; y < m_yBlockCnt; ++y) + for (size_t x = 0; x < m_xBlockCnt; ++x) writeBlock(x, y, si, srcNoData); } @@ -461,16 +475,16 @@ friend class Raster; } template - void writeBlock(int x, int y, SOURCE_ITER sourceBegin, + void writeBlock(size_t x, size_t y, SOURCE_ITER sourceBegin, ITER_VAL srcNoData) { - int xWidth = 0; + size_t xWidth = 0; if (x == m_xBlockCnt - 1) xWidth = m_xTotalSize % m_xBlockSize; if (xWidth == 0) xWidth = m_xBlockSize; - int yHeight = 0; + size_t yHeight = 0; if (y == m_yBlockCnt - 1) yHeight = m_yTotalSize % m_yBlockSize; if (yHeight == 0) @@ -480,11 +494,11 @@ friend class Raster; auto di = m_buf.begin(); // Go through rows copying data. Increment the destination iterator // by the width of the row. - for (int row = 0; row < yHeight; ++row) + for (size_t row = 0; row < yHeight; ++row) { // Find the offset location in the source container. - int wholeRowElts = m_xTotalSize * ((y * m_yBlockSize) + row); - int partialRowElts = m_xBlockSize * x; + size_t wholeRowElts = m_xTotalSize * ((y * m_yBlockSize) + row); + size_t partialRowElts = m_xBlockSize * x; auto si = sourceBegin + (wholeRowElts + partialRowElts); std::transform(si, si + xWidth, di, @@ -510,7 +524,10 @@ friend class Raster; // is valid, so we use m_xBlockSize instead of xWidth. di += m_xBlockSize; } - if (m_band->WriteBlock(x, y, m_buf.data()) != CPLE_None) + + // x and y are guaranteed to fit into an int + if (m_band->WriteBlock(static_cast(x), + static_cast(y), m_buf.data()) != CPLE_None) throw CantWriteBlock(); } @@ -520,9 +537,6 @@ friend class Raster; { m_band->GetStatistics(bApprox, bForce, minimum, maximum, mean, stddev); } - - - }; @@ -611,18 +625,20 @@ class PDAL_DLL Raster } catch (InvalidBand) { - std::stringstream oss; - oss << "Unable to get band " << nBand << " from raster '" << - m_filename << "'."; - m_errorMsg = oss.str(); + m_errorMsg = "Unable to get band " + std::to_string(nBand) + + " from raster '" + m_filename + "'."; return GDALError::InvalidBand; } + catch (BadBand) + { + m_errorMsg = "Unable to read band/block information from " + "raster '" + m_filename + "'."; + return GDALError::BadBand; + } catch (CantReadBlock) { - std::ostringstream oss; - oss << "Unable to read block for for raster '" << m_filename << - "'."; - m_errorMsg = oss.str(); + m_errorMsg = "Unable to read block for for raster '" + + m_filename + "'."; return GDALError::CantReadBlock; } return GDALError::None; @@ -644,58 +660,68 @@ class PDAL_DLL Raster { switch(m_bandType) { - case Dimension::Type::Unsigned8: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::Signed8: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::Unsigned16: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::Signed16: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::Unsigned32: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::Signed32: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::Unsigned64: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::Signed64: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::Float: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::Double: - Band(m_ds, nBand, m_dstNoData, name). - write(si, srcNoData); - break; - case Dimension::Type::None: - throw CantWriteBlock(); + case Dimension::Type::Unsigned8: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::Signed8: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::Unsigned16: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::Signed16: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::Unsigned32: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::Signed32: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::Unsigned64: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::Signed64: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::Float: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::Double: + Band(m_ds, nBand, m_dstNoData, name). + write(si, srcNoData); + break; + case Dimension::Type::None: + throw CantWriteBlock(); } } + catch (InvalidBand) + { + m_errorMsg = "Unable to get band " + std::to_string(nBand) + + " from raster '" + m_filename + "'."; + return GDALError::InvalidBand; + } + catch (BadBand) + { + m_errorMsg = "Unable to read band/block information from " + "raster '" + m_filename + "'."; + return GDALError::BadBand; + } catch (CantWriteBlock err) { - std::ostringstream oss; - oss << "Unable to write block for for raster '" << m_filename << - "'."; + m_errorMsg = "Unable to write block for for raster '" + + m_filename + "'."; if (err.what.size()) - oss << "\n" << err.what; - m_errorMsg = oss.str(); + m_errorMsg += "\n" + err.what; return GDALError::CantWriteBlock; } return GDALError::None; @@ -762,49 +788,37 @@ class PDAL_DLL Raster std::string const& filename() { return m_filename; } - void statistics(int nBand, - double* minimum, - double* maximum, - double* mean, - double* stddev, - int bApprox=TRUE, - int bForce=TRUE) const + void statistics(int nBand, double* minimum, double* maximum, double* mean, + double* stddev, int bApprox = TRUE, int bForce = TRUE) const { - Band(m_ds, nBand).statistics(minimum, maximum, mean, stddev, bApprox, bForce); + Band(m_ds, nBand).statistics(minimum, maximum, mean, stddev, + bApprox, bForce); } BOX2D bounds() const { - std::array coords; - pixelToCoord(height(), - width(), - coords); + pixelToCoord(height(), width(), coords); double maxx = coords[0]; double maxy = coords[1]; - pixelToCoord(0, - 0, - coords); + pixelToCoord(0, 0, coords); double minx = coords[0]; double miny = coords[1]; - BOX2D output(minx, miny, maxx, maxy); - return output; + return BOX2D(minx, miny, maxx, maxy); } BOX3D bounds(int nBand) const { - BOX2D box2 = bounds(); double minimum; double maximum; double mean; double stddev; statistics(nBand, &minimum, &maximum, &mean, &stddev); - BOX3D output(box2.minx, box2.miny, minimum, + return BOX3D(box2.minx, box2.miny, minimum, box2.maxx, box2.maxy, maximum); - return output; } private: diff --git a/pdal/KDIndex.hpp b/pdal/KDIndex.hpp index a3c1bee2f7..0aeb8d9ee6 100644 --- a/pdal/KDIndex.hpp +++ b/pdal/KDIndex.hpp @@ -101,25 +101,25 @@ class PDAL_DLL KD2Index : public KDIndex<2> throw pdal_error("KD2Index: point view missing 'Y' dimension."); } - PointId neighbor(double x, double y) + PointId neighbor(double x, double y) const { std::vector ids = neighbors(x, y, 1); return (ids.size() ? ids[0] : 0); } - PointId neighbor(PointId idx) + PointId neighbor(PointId idx) const { std::vector ids = neighbors(idx, 1); return (ids.size() ? ids[0] : 0); } - PointId neighbor(PointRef &point) + PointId neighbor(PointRef &point) const { std::vector ids = neighbors(point, 1); return (ids.size() ? ids[0] : 0); } - std::vector neighbors(double x, double y, point_count_t k) + std::vector neighbors(double x, double y, point_count_t k) const { k = (std::min)(m_buf.size(), k); std::vector output(k); @@ -135,7 +135,7 @@ class PDAL_DLL KD2Index : public KDIndex<2> return output; } - std::vector neighbors(PointId idx, point_count_t k) + std::vector neighbors(PointId idx, point_count_t k) const { double x = m_buf.getFieldAs(Dimension::Id::X, idx); double y = m_buf.getFieldAs(Dimension::Id::Y, idx); @@ -143,7 +143,7 @@ class PDAL_DLL KD2Index : public KDIndex<2> return neighbors(x, y, k); } - std::vector neighbors(PointRef &point, point_count_t k) + std::vector neighbors(PointRef &point, point_count_t k) const { double x = point.getFieldAs(Dimension::Id::X); double y = point.getFieldAs(Dimension::Id::Y); @@ -226,26 +226,26 @@ class PDAL_DLL KD3Index : public KDIndex<3> throw pdal_error("KD3Index: point view missing 'Z' dimension."); } - PointId neighbor(double x, double y, double z) + PointId neighbor(double x, double y, double z) const { std::vector ids = neighbors(x, y, z, 1); return (ids.size() ? ids[0] : 0); } - PointId neighbor(PointId idx) + PointId neighbor(PointId idx) const { std::vector ids = neighbors(idx, 1); return (ids.size() ? ids[0] : 0); } - PointId neighbor(PointRef &point) + PointId neighbor(PointRef &point) const { std::vector ids = neighbors(point, 1); return (ids.size() ? ids[0] : 0); } std::vector neighbors(double x, double y, double z, - point_count_t k) + point_count_t k) const { k = (std::min)(m_buf.size(), k); std::vector output(k); @@ -262,7 +262,7 @@ class PDAL_DLL KD3Index : public KDIndex<3> return output; } - std::vector neighbors(PointId idx, point_count_t k) + std::vector neighbors(PointId idx, point_count_t k) const { double x = m_buf.getFieldAs(Dimension::Id::X, idx); double y = m_buf.getFieldAs(Dimension::Id::Y, idx); @@ -271,7 +271,7 @@ class PDAL_DLL KD3Index : public KDIndex<3> return neighbors(x, y, z, k); } - std::vector neighbors(PointRef &point, point_count_t k) + std::vector neighbors(PointRef &point, point_count_t k) const { double x = point.getFieldAs(Dimension::Id::X); double y = point.getFieldAs(Dimension::Id::Y); diff --git a/pdal/PipelineWriter.cpp b/pdal/PipelineWriter.cpp index 91fb74be14..cd8ec2c369 100644 --- a/pdal/PipelineWriter.cpp +++ b/pdal/PipelineWriter.cpp @@ -62,10 +62,10 @@ std::string generateTag(Stage *stage, PipelineWriter::TagMap& tags) for (size_t i = 1; ; ++i) { tag = stage->getName() + std::to_string(i); + tag = Utils::replaceAll(tag, ".", "_"); if (!tagExists(tag)) break; } - tag = Utils::replaceAll(tag, ".", "_"); } return tag; } diff --git a/scripts/docker/ubuntu/Dockerfile b/scripts/docker/ubuntu/Dockerfile index 9b3e39b39a..a302a6087e 100644 --- a/scripts/docker/ubuntu/Dockerfile +++ b/scripts/docker/ubuntu/Dockerfile @@ -73,7 +73,7 @@ RUN apt-get update && DEBIAN_FRONTEND=noninteractive apt-get install -y --fix-m RUN git clone https://github.com/LASzip/LASzip.git laszip; \ cd laszip; \ - git checkout 3.1.1; \ + git checkout 3.4.1; \ cmake \ -G Ninja \ -DCMAKE_INSTALL_PREFIX=/usr/ \ @@ -292,10 +292,10 @@ COPY --from=builder /build/usr/lib/ /usr/lib/ COPY --from=builder /build/usr/include/ /usr/include/ RUN \ - curl -LOs http://download.osgeo.org/proj/proj-datumgrid-1.8.zip && unzip -j -u -o proj-datumgrid-1.8.zip -d /usr/share/proj; \ - curl -LOs http://download.osgeo.org/proj/proj-datumgrid-europe-1.2.zip && unzip -j -u -o proj-datumgrid-europe-1.2.zip -d /usr/share/proj; \ - curl -LOs http://download.osgeo.org/proj/proj-datumgrid-oceania-1.0.zip && unzip -j -u -o proj-datumgrid-oceania-1.0.zip -d /usr/share/proj; \ - curl -LOs http://download.osgeo.org/proj/proj-datumgrid-world-1.0.zip && unzip -j -u -o proj-datumgrid-world-1.0.zip -d /usr/share/proj; \ - curl -LOs http://download.osgeo.org/proj/proj-datumgrid-north-america-1.2.zip && unzip -j -u -o proj-datumgrid-north-america-1.2.zip -d /usr/share/proj; + curl -LOs http://download.osgeo.org/proj/proj-datumgrid-1.8.zip && unzip -j -u -o proj-datumgrid-1.8.zip -d /usr/share/proj && rm proj-datumgrid-1.8.zip ; \ + curl -LOs http://download.osgeo.org/proj/proj-datumgrid-europe-1.2.zip && unzip -j -u -o proj-datumgrid-europe-1.2.zip -d /usr/share/proj && rm proj-datumgrid-europe-1.2.zip ; \ + curl -LOs http://download.osgeo.org/proj/proj-datumgrid-oceania-1.0.zip && unzip -j -u -o proj-datumgrid-oceania-1.0.zip -d /usr/share/proj && rm proj-datumgrid-oceania-1.0.zip; \ + curl -LOs http://download.osgeo.org/proj/proj-datumgrid-world-1.0.zip && unzip -j -u -o proj-datumgrid-world-1.0.zip -d /usr/share/proj && rm proj-datumgrid-world-1.0.zip; \ + curl -LOs http://download.osgeo.org/proj/proj-datumgrid-north-america-1.2.zip && unzip -j -u -o proj-datumgrid-north-america-1.2.zip -d /usr/share/proj && rm proj-datumgrid-north-america-1.2.zip; diff --git a/test/data/ply/issue_2421.ply b/test/data/ply/issue_2421.ply index 047f8b87a4..92145c32cf 100644 --- a/test/data/ply/issue_2421.ply +++ b/test/data/ply/issue_2421.ply @@ -7,4 +7,4 @@ property float64 y property float64 z property int32 i end_header -1.23457 12345.67890 1234567890.12345 12345 +1.23457 12345.67890 1234567890.12345 1234567890 diff --git a/test/data/ply/sized_dims.ply b/test/data/ply/sized_dims.ply new file mode 100644 index 0000000000..bcfbb9281e --- /dev/null +++ b/test/data/ply/sized_dims.ply @@ -0,0 +1,10 @@ +ply +format ascii 1.0 +comment Generated by PDAL +element vertex 1 +property int8 x +property uint16 y +property int32 j +property float64 i +end_header +1 12346 12345 1234567890.00000 diff --git a/test/data/ply/unsized_dims.ply b/test/data/ply/unsized_dims.ply new file mode 100644 index 0000000000..04b4055834 --- /dev/null +++ b/test/data/ply/unsized_dims.ply @@ -0,0 +1,10 @@ +ply +format ascii 1.0 +comment Generated by PDAL +element vertex 1 +property ushort y +property char x +property double i +property int j +end_header +12346 1 1234567890.00000 12345 diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index a6bb8ae355..cb0d26f47c 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -52,6 +52,12 @@ PDAL_ADD_TEST(pdal_options_test ) PDAL_ADD_TEST(pdal_pipeline_manager_test FILES PipelineManagerTest.cpp) +PDAL_ADD_TEST(pdal_pipeline_writer_test + FILES + PipelineWriterTest.cpp + INCLUDES + ${PDAL_JSONCPP_INCLUDE_DIR} +) PDAL_ADD_TEST(pdal_plugin_manager_test FILES PluginManagerTest.cpp) PDAL_ADD_TEST(pdal_point_view_test FILES PointViewTest.cpp) PDAL_ADD_TEST(pdal_point_table_test FILES PointTableTest.cpp) @@ -221,6 +227,7 @@ PDAL_ADD_TEST(pdal_filters_crop_test PDAL_ADD_TEST(pdal_filters_decimation_test FILES filters/DecimationFilterTest.cpp) PDAL_ADD_TEST(pdal_filters_delaunay_test FILES filters/DelaunayFilterTest.cpp) +PDAL_ADD_TEST(pdal_filters_covariancefeatures_test FILES filters/CovarianceFeaturesTest.cpp) PDAL_ADD_TEST(pdal_filters_divider_test FILES filters/DividerFilterTest.cpp) PDAL_ADD_TEST(pdal_filters_mongoexpression_test FILES diff --git a/test/unit/PipelineWriterTest.cpp b/test/unit/PipelineWriterTest.cpp new file mode 100644 index 0000000000..c54d2cda4b --- /dev/null +++ b/test/unit/PipelineWriterTest.cpp @@ -0,0 +1,69 @@ +/****************************************************************************** +* Copyright (c) 2019, Michael P. Gerlek (mpg@flaxen.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 + +#include + +#include "Support.hpp" + +#include +#include +#include + +using namespace pdal; + +// Make sure we handle duplicate stages properly. +TEST(PipelineManagerTest, issue_2458) +{ + std::string in = R"( + [ + "in.las", + "in2.las", + "out.las" + ] + )"; + + PipelineManager mgr; + std::istringstream iss(in); + mgr.readPipeline(iss); + + std::ostringstream oss; + PipelineWriter::writePipeline(mgr.getStage(), oss); + + std::string out = oss.str(); + EXPECT_TRUE(out.find("readers_las1") != std::string::npos); + EXPECT_TRUE(out.find("readers_las2") != std::string::npos); + EXPECT_TRUE(out.find("writers_las1") != std::string::npos); +} diff --git a/test/unit/filters/CovarianceFeaturesTest.cpp b/test/unit/filters/CovarianceFeaturesTest.cpp new file mode 100644 index 0000000000..fae3f94888 --- /dev/null +++ b/test/unit/filters/CovarianceFeaturesTest.cpp @@ -0,0 +1,184 @@ +/****************************************************************************** +* Copyright (c) 2019, Helix Re Inc. nicolas@helix.re +* +* 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 Helix Re Inc. 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 + +#include +#include +#include +#include + +#include "Support.hpp" + +namespace pdal { + +TEST(DimensionalityTest, Linearity) +{ + using namespace Dimension; + + PointTable table; + table.layout()->registerDims({Id::X, Id::Y, Id::Z}); + + BufferReader bufferReader; + CovarianceFeaturesFilter filter; + Options ops; + ops.add("knn", 3); + filter.setInput(bufferReader); + filter.setOptions(ops); + filter.prepare(table); + + PointViewPtr view(new PointView(table)); + view->setField(Id::X, 0,0); + view->setField(Id::X, 1, 0); + view->setField(Id::X, 2, 0); + view->setField(Id::Y, 0, 0); + view->setField(Id::Y, 1, 0); + view->setField(Id::Y, 2, 0); + view->setField(Id::Z, 0, 0); + view->setField(Id::Z, 1, 1); + view->setField(Id::Z, 2, 2); + bufferReader.addView(view); + + PointViewSet viewSet = filter.execute(table); + PointViewPtr outView = *viewSet.begin(); + + Dimension::Id planarity = table.layout()->findDim("Planarity"); + Dimension::Id verticality = table.layout()->findDim("Verticality"); + Dimension::Id linearity = table.layout()->findDim("Linearity"); + Dimension::Id scattering = table.layout()->findDim("Scattering"); + + for (point_count_t i =0; i < outView->size(); i++) + { + ASSERT_FLOAT_EQ(outView->getFieldAs(linearity, i) ,1); + ASSERT_FLOAT_EQ(outView->getFieldAs(verticality, i) ,1); + ASSERT_FLOAT_EQ(outView->getFieldAs(planarity, i) ,0); + ASSERT_FLOAT_EQ(outView->getFieldAs(scattering, i) ,0); + } +} + +TEST(DimensionalityTest, Planarity) +{ + using namespace Dimension; + + PointTable table; + table.layout()->registerDims({Id::X, Id::Y, Id::Z}); + + BufferReader bufferReader; + CovarianceFeaturesFilter filter; + Options ops; + filter.setInput(bufferReader); + filter.prepare(table); + + PointViewPtr view(new PointView(table)); + view->setField(Id::X, 0,0); + view->setField(Id::X, 1, 1); + view->setField(Id::X, 2, 1); + view->setField(Id::X, 3, 0); + view->setField(Id::Y, 0, 0); + view->setField(Id::Y, 1, 0); + view->setField(Id::Y, 2, 1); + view->setField(Id::Y, 3, 1); + view->setField(Id::Z, 0, 0); + view->setField(Id::Z, 1, 0); + view->setField(Id::Z, 2, 0); + view->setField(Id::Z, 3, 0); + + bufferReader.addView(view); + + + PointViewSet viewSet = filter.execute(table); + PointViewPtr outView = *viewSet.begin(); + + Dimension::Id planarity = table.layout()->findDim("Planarity"); + Dimension::Id verticality = table.layout()->findDim("Verticality"); + Dimension::Id linearity = table.layout()->findDim("Linearity"); + Dimension::Id scattering = table.layout()->findDim("Scattering"); + + for (point_count_t i =0; i < outView->size(); i++) + { + ASSERT_LE(outView->getFieldAs(linearity, i) ,0.5); + ASSERT_FLOAT_EQ(outView->getFieldAs(verticality, i) ,0); + ASSERT_FLOAT_EQ(outView->getFieldAs(planarity, i) ,1); + ASSERT_FLOAT_EQ(outView->getFieldAs(scattering, i) ,0); + } +} + +TEST(DimensionalityTest, Scattering) +{ + using namespace Dimension; + + PointTable table; + table.layout()->registerDims({Id::X, Id::Y, Id::Z}); + + BufferReader bufferReader; + CovarianceFeaturesFilter filter; + Options ops; + filter.setInput(bufferReader); + filter.prepare(table); + + PointViewPtr view(new PointView(table)); + view->setField(Id::X, 0,0); + view->setField(Id::X, 1, 1); + view->setField(Id::X, 2, 1); + view->setField(Id::X, 3, 0); + view->setField(Id::Y, 0, 0); + view->setField(Id::Y, 1, 0); + view->setField(Id::Y, 2, 1); + view->setField(Id::Y, 3, 1); + view->setField(Id::Z, 0, 0); + view->setField(Id::Z, 1, 0); + view->setField(Id::Z, 2, 1); + view->setField(Id::Z, 3, 2); + + bufferReader.addView(view); + + + PointViewSet viewSet = filter.execute(table); + PointViewPtr outView = *viewSet.begin(); + + Dimension::Id planarity = table.layout()->findDim("Planarity"); + Dimension::Id verticality = table.layout()->findDim("Verticality"); + Dimension::Id linearity = table.layout()->findDim("Linearity"); + Dimension::Id scattering = table.layout()->findDim("Scattering"); + + for (point_count_t i =0; i < outView->size(); i++) + { + ASSERT_LE(outView->getFieldAs(linearity, i) ,0.5); + ASSERT_GE(outView->getFieldAs(verticality, i) ,0.5); + ASSERT_LE(outView->getFieldAs(planarity, i) ,1); + ASSERT_GE(outView->getFieldAs(scattering, i) ,0.1); + } +} + +} \ No newline at end of file diff --git a/test/unit/io/PlyWriterTest.cpp b/test/unit/io/PlyWriterTest.cpp index 3b96b4d43b..d53024ccb2 100644 --- a/test/unit/io/PlyWriterTest.cpp +++ b/test/unit/io/PlyWriterTest.cpp @@ -138,6 +138,9 @@ TEST(PlyWriter, mesh) // Make sure that you can't use precision with ascii output. TEST(PlyWriter, precisionException) { + std::string outfile(Support::temppath("out.ply")); + FileUtils::deleteFile(outfile); + Options readerOptions; readerOptions.add("count", 750); readerOptions.add("mode", "random"); @@ -145,7 +148,7 @@ TEST(PlyWriter, precisionException) reader.setOptions(readerOptions); Options writerOptions; - writerOptions.add("filename", Support::temppath("out.ply")); + writerOptions.add("filename", outfile); writerOptions.add("storage_mode", "little endian"); writerOptions.add("precision", 3); PlyWriter writer; @@ -175,7 +178,7 @@ TEST(PlyWriter, issue_2421) v->setField(Dimension::Id::X, 0, 1.23456789012345); v->setField(Dimension::Id::Y, 0, 12345.6789012345); v->setField(Dimension::Id::Z, 0, 1234567890.12345); - v->setField(iid, 0, 12345); + v->setField(iid, 0, 1234567890); BufferReader r; r.addView(v); @@ -194,4 +197,64 @@ TEST(PlyWriter, issue_2421) EXPECT_TRUE(Support::compare_text_files(outfile, referenceFile)); } +TEST(PlyWriter, dimtypes) +{ + std::string outfile(Support::temppath("out.ply")); + + PointTable t; + + t.layout()->registerDim(Dimension::Id::X); + t.layout()->registerDim(Dimension::Id::Y); + t.layout()->registerDim(Dimension::Id::Z); + Dimension::Id iid = t.layout()->assignDim("I", Dimension::Type::Double); + Dimension::Id jid = t.layout()->assignDim("J", Dimension::Type::Double); + + PointViewPtr v(new PointView(t)); + v->setField(Dimension::Id::X, 0, 1.23456789012345); + v->setField(Dimension::Id::Y, 0, 12345.6789012345); + v->setField(Dimension::Id::Z, 0, 1234567890.12345); + v->setField(iid, 0, 1234567890); + v->setField(jid, 0, 12345); + + BufferReader r; + r.addView(v); + + { + PlyWriter w; + Options wo; + wo.add("filename", outfile); + wo.add("storage_mode", "ascii"); + wo.add("precision", 5); + wo.add("dims", "X=int8_t,Y=ushort,J=int,I=double"); + w.setInput(r); + w.setOptions(wo); + + FileUtils::deleteFile(outfile); + w.prepare(t); + w.execute(t); + + std::string referenceFile(Support::datapath("ply/sized_dims.ply")); + EXPECT_TRUE(Support::compare_text_files(outfile, referenceFile)); + } + + { + PlyWriter w; + Options wo; + wo.add("filename", outfile); + wo.add("storage_mode", "ascii"); + wo.add("precision", 5); + wo.add("sized_types", false); + wo.add("dims", "Y=ushort,X=int8_t,I,J=int"); + w.setInput(r); + w.setOptions(wo); + + FileUtils::deleteFile(outfile); + w.prepare(t); + w.execute(t); + + std::string referenceFile(Support::datapath("ply/unsized_dims.ply")); + EXPECT_TRUE(Support::compare_text_files(outfile, referenceFile)); + } +} + } // namespace pdal