From 8b48f379c442409c6462f14aaed16a46644266eb Mon Sep 17 00:00:00 2001 From: Howard Butler Date: Wed, 12 Apr 2017 12:48:30 -0500 Subject: [PATCH 01/57] arbiter refresh --- vendor/arbiter/arbiter.cpp | 37 ++++++++++--- vendor/arbiter/arbiter.hpp | 108 ++++++++++++++++++++++++++++--------- 2 files changed, 114 insertions(+), 31 deletions(-) diff --git a/vendor/arbiter/arbiter.cpp b/vendor/arbiter/arbiter.cpp index 349878c875..f39bdf5a71 100644 --- a/vendor/arbiter/arbiter.cpp +++ b/vendor/arbiter/arbiter.cpp @@ -79,7 +79,27 @@ namespace const std::size_t httpRetryCount(8); #endif - Json::Value getConfig() + // Merge B into A, without overwriting any keys from A. + Json::Value merge(const Json::Value& a, const Json::Value& b) + { + Json::Value out(a); + + if (!b.isNull()) + { + for (const auto& key : b.getMemberNames()) + { + // If A doesn't have this key, then set it to B's value. + // If A has the key but it's an object, then recursively merge. + // Otherwise A already has a value here that we won't overwrite. + if (!out.isMember(key)) out[key] = b[key]; + else if (out[key].isObject()) merge(out[key], b[key]); + } + } + + return out; + } + + Json::Value getConfig(const Json::Value& in) { Json::Value config; std::string path("~/.arbiter/config.json"); @@ -93,20 +113,22 @@ namespace ss >> config; } - return config; + return merge(in, config); } } -Arbiter::Arbiter() : Arbiter(getConfig()) { } +Arbiter::Arbiter() : Arbiter(Json::nullValue) { } -Arbiter::Arbiter(const Json::Value& json) +Arbiter::Arbiter(const Json::Value& in) : m_drivers() #ifdef ARBITER_CURL - , m_pool(new http::Pool(concurrentHttpReqs, httpRetryCount, json)) + , m_pool(new http::Pool(concurrentHttpReqs, httpRetryCount, getConfig(in))) #endif { using namespace drivers; + const Json::Value json(getConfig(in)); + auto fs(Fs::create(json["file"])); if (fs) m_drivers[fs->type()] = std::move(fs); @@ -1747,7 +1769,7 @@ S3::AuthFields S3::Auth::fields() const m_access = creds["AccessKeyId"].asString(); m_hidden = creds["SecretAccessKey"].asString(); m_token = creds["Token"].asString(); - m_expiration.reset(new Time(creds["Expiration"].asString())); + m_expiration.reset(new Time(creds["Expiration"].asString(), arbiter::Time::iso8601)); if (*m_expiration - now < reauthSeconds) { @@ -2619,7 +2641,7 @@ namespace PutData* in) { const std::size_t fullBytes( - std::min( + (std::min)( size * num, in->data.size() - in->offset)); std::memcpy(out, in->data.data() + in->offset, fullBytes); @@ -4022,6 +4044,7 @@ int64_t Time::operator-(const Time& other) const #endif #include +#include #ifdef ARBITER_CUSTOM_NAMESPACE namespace ARBITER_CUSTOM_NAMESPACE diff --git a/vendor/arbiter/arbiter.hpp b/vendor/arbiter/arbiter.hpp index bb00983ec6..6c28dd886f 100644 --- a/vendor/arbiter/arbiter.hpp +++ b/vendor/arbiter/arbiter.hpp @@ -1,7 +1,7 @@ /// Arbiter amalgamated header (https://github.com/connormanning/arbiter). /// It is intended to be used with #include "arbiter.hpp" -// Git SHA: 5717b5f7b450db996b4a44d4865671f925e2ff50 +// Git SHA: 1e9cadf66ddc1bc5ea7d4e044136f6072b8b305a // ////////////////////////////////////////////////////////////////////// // Beginning of content of file: LICENSE @@ -48,6 +48,43 @@ SOFTWARE. #define ARBITER_CUSTOM_NAMESPACE pdal #define ARBITER_EXTERNAL_JSON +// ////////////////////////////////////////////////////////////////////// +// Beginning of content of file: arbiter/util/exports.hpp +// ////////////////////////////////////////////////////////////////////// + +#pragma once + +#if defined(_WIN32) || defined(WIN32) || defined(_MSC_VER) +#define ARBITER_WINDOWS +#endif + +#ifndef ARBITER_DLL +#if defined(ARBITER_WINDOWS) +#if defined(ARBITER_DLL_EXPORT) +# define ARBITER_DLL __declspec(dllexport) +#elif defined(PDAL_DLL_IMPORT) +# define ARBITER_DLL __declspec(dllimport) +#else +# define ARBITER_DLL +#endif +#else +# if defined(USE_GCC_VISIBILITY_FLAG) +# define ARBITER_DLL __attribute__ ((visibility("default"))) +# else +# define ARBITER_DLL +# endif +#endif +#endif + +// ////////////////////////////////////////////////////////////////////// +// End of content of file: arbiter/util/exports.hpp +// ////////////////////////////////////////////////////////////////////// + + + + + + // ////////////////////////////////////////////////////////////////////// // Beginning of content of file: arbiter/util/types.hpp // ////////////////////////////////////////////////////////////////////// @@ -156,6 +193,8 @@ class Response #ifndef ARBITER_IS_AMALGAMATION #include +#include + #ifndef ARBITER_EXTERNAL_JSON #include @@ -169,7 +208,11 @@ class Response #include #endif -class curl_slist; +#ifdef ARBITER_CURL +#include +#endif + +struct curl_slist; #ifdef ARBITER_CUSTOM_NAMESPACE namespace ARBITER_CUSTOM_NAMESPACE @@ -185,7 +228,7 @@ namespace http class Pool; -class Curl +class ARBITER_DLL Curl { friend class Pool; @@ -276,6 +319,7 @@ class Curl #include #include +#include #ifndef ARBITER_EXTERNAL_JSON #include @@ -312,9 +356,9 @@ std::string buildQueryString(const http::Query& query); /** @cond arbiter_internal */ -class Pool; +class ARBITER_DLL Pool; -class Resource +class ARBITER_DLL Resource { public: Resource(Pool& pool, Curl& curl, std::size_t id, std::size_t retry); @@ -352,7 +396,7 @@ class Resource http::Response exec(std::function f); }; -class Pool +class ARBITER_DLL Pool { // Only HttpResource may release. friend class Resource; @@ -451,6 +495,12 @@ Contents parse(const std::string& s); #include #include +#ifndef ARBITER_IS_AMALGAMATION +#include +#endif + + + #ifdef ARBITER_CUSTOM_NAMESPACE namespace ARBITER_CUSTOM_NAMESPACE { @@ -459,7 +509,7 @@ namespace ARBITER_CUSTOM_NAMESPACE namespace arbiter { -class Time +class ARBITER_DLL Time { public: static const std::string iso8601; @@ -467,9 +517,9 @@ class Time static const std::string dateNoSeparators; Time(); - Time(const std::string& s, const std::string& format = iso8601); + Time(const std::string& s, const std::string& format = "%Y-%m-%dT%H:%M:%SZ"); - std::string str(const std::string& format = iso8601) const; + std::string str(const std::string& format = "%Y-%m-%dT%H:%M:%SZ") const; // Return value is in seconds. int64_t operator-(const Time& other) const; @@ -654,7 +704,9 @@ typedef std::map> DriverMap; #endif - +#ifndef ARBITER_IS_AMALGAMATION +#include +#endif #ifdef ARBITER_EXTERNAL_JSON #include @@ -679,20 +731,20 @@ class Arbiter; namespace fs { /** @brief Returns true if created, false if already existed. */ - bool mkdirp(std::string dir); + ARBITER_DLL bool mkdirp(std::string dir); /** @brief Returns true if removed, otherwise false. */ - bool remove(std::string filename); + ARBITER_DLL bool remove(std::string filename); /** @brief Performs tilde expansion to a fully-qualified path, if possible. */ - std::string expandTilde(std::string path); + ARBITER_DLL std::string expandTilde(std::string path); /** @brief Get temporary path from environment. */ - std::string getTempPath(); + ARBITER_DLL std::string getTempPath(); /** @brief Resolve a possible wildcard path. */ - std::vector glob(std::string path); + ARBITER_DLL std::vector glob(std::string path); /** @brief A scoped local filehandle for a possibly remote path. * @@ -702,7 +754,7 @@ namespace fs * * See Arbiter::getLocalHandle for details about construction. */ - class LocalHandle + class ARBITER_DLL LocalHandle { friend class arbiter::Arbiter; @@ -746,9 +798,11 @@ namespace drivers { /** @brief Local filesystem driver. */ -class Fs : public Driver +class ARBITER_DLL Fs : public Driver { public: + Fs() { } + static std::unique_ptr create(const Json::Value& json); virtual std::string type() const override { return "file"; } @@ -3824,6 +3878,12 @@ std::string encodeAsHex(const std::string& data); #include #include #include +#include + +#ifndef ARBITER_IS_AMALGAMATION +#include +#endif + #ifdef ARBITER_CUSTOM_NAMESPACE namespace ARBITER_CUSTOM_NAMESPACE @@ -3839,7 +3899,7 @@ namespace util /** Returns @p path, less any trailing glob indicators (one or two * asterisks) as well as any possible trailing slash. */ - std::string stripPostfixing(std::string path); + ARBITER_DLL std::string stripPostfixing(std::string path); /** Returns the portion of @p fullPath following the last instance of the * character `/`, if any instances exist aside from possibly the delimiter @@ -3851,20 +3911,20 @@ namespace util * logic above, thus the innermost directory in the full path will be * returned. */ - std::string getBasename(std::string fullPath); + ARBITER_DLL std::string getBasename(std::string fullPath); /** Returns everything besides the basename, as determined by `getBasename`. * For file paths, this corresponds to the directory path above the file. * For directory paths, this corresponds to all directories above the * innermost directory. */ - std::string getNonBasename(std::string fullPath); + ARBITER_DLL std::string getNonBasename(std::string fullPath); /** @cond arbiter_internal */ - inline bool isSlash(char c) { return c == '/' || c == '\\'; } + ARBITER_DLL inline bool isSlash(char c) { return c == '/' || c == '\\'; } /** Returns true if the last character is an asterisk. */ - inline bool isGlob(std::string path) + ARBITER_DLL inline bool isGlob(std::string path) { return path.size() && path.back() == '*'; } @@ -4631,7 +4691,7 @@ class Endpoint #endif #ifndef ARBITER_IS_AMALGAMATION - +#include #include #include #include @@ -4675,7 +4735,7 @@ namespace http { class Pool; } * * All Arbiter operations are thread-safe except unless otherwise noted. */ -class Arbiter +class ARBITER_DLL Arbiter { public: /** Construct a basic Arbiter with only drivers the don't require From 4f6babc386cad5956ffeb58a3934e20a6d4277cd Mon Sep 17 00:00:00 2001 From: Howard Butler Date: Wed, 12 Apr 2017 12:56:43 -0500 Subject: [PATCH 02/57] another arbiter refresh --- vendor/arbiter/arbiter.cpp | 6 +++--- vendor/arbiter/arbiter.hpp | 7 ++++++- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/vendor/arbiter/arbiter.cpp b/vendor/arbiter/arbiter.cpp index f39bdf5a71..797358fdca 100644 --- a/vendor/arbiter/arbiter.cpp +++ b/vendor/arbiter/arbiter.cpp @@ -1250,7 +1250,7 @@ std::unique_ptr Http::tryGetSize(std::string path) const std::unique_ptr size; auto http(m_pool.acquire()); - Response res(http.head(path)); + Response res(http.head(typedPath(path))); if (res.ok() && res.headers().count("Content-Length")) { @@ -1322,7 +1322,7 @@ bool Http::get( bool good(false); auto http(m_pool.acquire()); - Response res(http.get(path, headers, query)); + Response res(http.get(typedPath(path), headers, query)); if (res.ok()) { @@ -1341,7 +1341,7 @@ void Http::put( { auto http(m_pool.acquire()); - if (!http.put(path, data, headers, query).ok()) + if (!http.put(typedPath(path), data, headers, query).ok()) { throw ArbiterError("Couldn't HTTP PUT to " + path); } diff --git a/vendor/arbiter/arbiter.hpp b/vendor/arbiter/arbiter.hpp index 6c28dd886f..50bf7a6dc1 100644 --- a/vendor/arbiter/arbiter.hpp +++ b/vendor/arbiter/arbiter.hpp @@ -1,7 +1,7 @@ /// Arbiter amalgamated header (https://github.com/connormanning/arbiter). /// It is intended to be used with #include "arbiter.hpp" -// Git SHA: 1e9cadf66ddc1bc5ea7d4e044136f6072b8b305a +// Git SHA: 311b81e3dcb5e8c55b2148f7eba46a12cff66914 // ////////////////////////////////////////////////////////////////////// // Beginning of content of file: LICENSE @@ -994,6 +994,11 @@ class Http : public Driver return get(path, data, http::Headers(), http::Query()); } + std::string typedPath(const std::string& p) const + { + return type() + "://" + p; + } + http::Pool& m_pool; }; From 2abc9955f690d781c3b1face9d30f5373d0464c1 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Wed, 12 Apr 2017 14:46:34 -0500 Subject: [PATCH 03/57] Initial check-in for safety. --- filters/GreedyProjection.cpp | 1870 ++++++++++++++++++++++++++++++++++ filters/GreedyProjection.hpp | 459 +++++++++ pdal/Dimension.json | 15 + 3 files changed, 2344 insertions(+) create mode 100644 filters/GreedyProjection.cpp create mode 100644 filters/GreedyProjection.hpp diff --git a/filters/GreedyProjection.cpp b/filters/GreedyProjection.cpp new file mode 100644 index 0000000000..1e423ebbdb --- /dev/null +++ b/filters/GreedyProjection.cpp @@ -0,0 +1,1870 @@ +/* + * Software License Agreement (BSD License) + * + * Copyright (c) 2010, Willow Garage, Inc. + * 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 Willow Garage, 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. + * + * $Id$ + * + */ + +#include + +#include + +#include "GreedyProjection.hpp" + +namespace pdal +{ + +void GreedyProjection::addArgs(ProgramArgs& args) +{ + args.add("multiplier", "Nearest neighbor distance multiplier", + mu_).setPositional(); + args.add("radius", "Search radius for neighbors", + search_radius_).setPositional(); + args.add("num_neighbors", "Number of nearest neighbors to consider", + nnn_, 100); + args.add("min_angle", "Minimum angle for created triangles", + minimum_angle_, M_PI / 18); // 10 degrees default + args.add("max_angle", "Maximum angle for created triangles", + maximum_angle_, 2 * M_PI / 3); // 120 degrees default + args.add("eps_angle", "Max normal difference angle for triangulation " + "consideration", eps_angle_, M_PI / 4); +} + +void GreedyProjection::initialize() +{ + if (search_radius_ <= 0) + throwError("Invalid search radius of '" + + std::to_string(search_radius_) + "'. Must be greater than 0."); + if (mu_ <= 0) + throwError("Invalid distance multiplier of '" + + std::to_string(mu_) + "'. Must be greater than 0."); +} + +Eigen::Vector3d GreedyProjection::getCoord(PointId id) +{ + assert(view_); + + return Eigen::Vector3d( + view_->getFieldAs(Dimension::Id::X, id), + view_->getFieldAs(Dimension::Id::Y, id), + view_->getFieldAs(Dimension::Id::Z, id) + ); +} + + +Eigen::Vector3d GreedyProjection::getNormalCoord(PointId id) +{ + assert(view_); + + return Eigen::Vector3d( + view_->getFieldAs(Dimension::Id::NormalX, id), + view_->getFieldAs(Dimension::Id::NormalY, id), + view_->getFieldAs(Dimension::Id::NormalZ, id) + ); +} + + +void GreedyProjection::addTriangle(PointId a, PointId b, PointId c) +{ + grid_.emplace_back(a, b, c); +} + + +void GreedyProjection::filter(PointView& view) +{ + KD3Index tree(view); + + view_ = &view; + const double sqr_mu = mu_ * mu_; + const double sqr_max_edge = search_radius_*search_radius_; + + nnn_ = (int)std::min((point_count_t)nnn_, view.size()); + + // Variables to hold the results of nearest neighbor searches + std::vector nnIdx(nnn_); + std::vector sqrDists(nnn_); + + // current number of connected components + int part_index = 0; + + // 2D coordinates of points + const Eigen::Vector2f uvn_nn_qp_zero = Eigen::Vector2f::Zero(); + Eigen::Vector2f uvn_current; + Eigen::Vector2f uvn_prev; + Eigen::Vector2f uvn_next; + + // initializing fields + already_connected_ = false; // see declaration for comments :P + + // initializing states and fringe neighbors + part_.clear (); + state_.clear (); + source_.clear (); + ffn_.clear (); + sfn_.clear (); + part_.resize(view.size (), -1); // indices of point's part + state_.resize(view.size (), GP3Type::FREE); + source_.resize(view.size (), -1); + ffn_.resize(view.size (), -1); + sfn_.resize(view.size (), -1); + fringe_queue_.clear (); + int fqIdx = 0; // current fringe's index in the queue to be processed + + // Avoiding NaN coordinates if needed + /** + if (!input_->is_dense) + { + // Skip invalid points from the indices list + for (std::vector::const_iterator it = indices_->begin ();o + it != indices_->end (); ++it) + if (!pcl_isfinite (input_->points[*it].x) || + !pcl_isfinite (input_->points[*it].y) || + !pcl_isfinite (input_->points[*it].z)) + state_[*it] = GP3Type::NONE; + } + **/ + + // Saving coordinates and point to index mapping + /** + coords_.clear (); + coords_.reserve (indices_->size ()); + 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()); + point2index[(*indices_)[cp]] = cp; + } + **/ + + // Initializing + PointId isFree = 0; + bool done = false; + int is_free=0, nr_parts=0, increase_nnn4fn=0, increase_nnn4s=0, increase_dist=0, nr_touched = 0; + bool is_fringe; + angles_.resize(nnn_); + std::vector > uvn_nn (nnn_); + Eigen::Vector2f uvn_s; + + // iterating through fringe points and finishing them until everything is done + while (!done) + { + R_ = is_free; + if (state_[R_] == GP3Type::FREE) + { + state_[R_] = GP3Type::NONE; +//ABELL - What is part_index? + part_[R_] = part_index++; + + // creating starting triangle + tree.knnSearch(R_, nnn_, &nnIdx, &sqrDists); + double sqr_dist_threshold = + (std::min)(sqr_max_edge, sqr_mu * sqrDists[1]); + +/** + // Search tree returns indices into the original cloud, but we are working with indices. TODO: make that optional! + for (int i = 1; i < nnn_; i++) + { + //if (point2index[nnIdx[i]] == -1) + // std::cerr << R_ << " [" << indices_->at (R_) << "] " << i << ": " << nnIdx[i] << " / " << point2index[nnIdx[i]] << std::endl; + nnIdx[i] = point2index[nnIdx[i]]; + } +**/ + + // Get the normal estimate at the current point + const Eigen::Vector3d nc(getNormalCoord(R_)); + + // Get a coordinate system that lies on a plane defined by its normal + v_ = nc.unitOrthogonal (); + u_ = nc.cross (v_); + + // Projecting point onto the surface + const Eigen::Vector3d coord(getCoord(R_)); + float dist = nc.dot (coord); + proj_qp_ = coord - dist * nc; + + // Converting coords, calculating angles and saving the + // projected near boundary edges + int nr_edge = 0; + std::vector doubleEdges; + // nearest neighbor with index 0 is the query point R_ itself + //ABELL - Need to check that the above is true. + for (int i = 1; i < nnn_; i++) + { + // Transforming coordinates + tmp_ = getCoord(nnIdx[i]) - proj_qp_; + uvn_nn[i][0] = tmp_.dot(u_); + uvn_nn[i][1] = tmp_.dot(v_); + // Computing the angle between each neighboring point and + // the query point itself + angles_[i].angle = atan2(uvn_nn[i][1], uvn_nn[i][0]); + // initializing angle descriptors + angles_[i].index = nnIdx[i]; + if ( + (state_[nnIdx[i]] == GP3Type::COMPLETED) || + (state_[nnIdx[i]] == GP3Type::BOUNDARY) || + (state_[nnIdx[i]] == GP3Type::NONE) || + (nnIdx[i] == -1) /// NOTE: discarding NaN points and those that are not in indices_ + || (sqrDists[i] > sqr_dist_threshold) + ) + angles_[i].visible = false; + else + angles_[i].visible = true; + // Saving the edges between nearby boundary points + if ((state_[nnIdx[i]] == GP3Type::FRINGE) || + (state_[nnIdx[i]] == GP3Type::BOUNDARY)) + { + doubleEdge e; + e.index = i; + nr_edge++; + tmp_ = getCoord(ffn_[nnIdx[i]]) - proj_qp_; + e.first[0] = tmp_.dot(u_); + e.first[1] = tmp_.dot(v_); + tmp_ = getCoord(sfn_[nnIdx[i]]) - proj_qp_; + e.second[0] = tmp_.dot(u_); + e.second[1] = tmp_.dot(v_); + doubleEdges.push_back(e); + } + } + angles_[0].visible = false; + + // Verify the visibility of each potential new vertex + // nearest neighbor with index 0 is the query point R_ itself + // ABELL - Check the above. + for (int i = 1; i < nnn_; i++) + if ((angles_[i].visible) && (ffn_[R_] != nnIdx[i]) && + (sfn_[R_] != nnIdx[i])) + { + bool visibility = true; + for (int j = 0; j < nr_edge; j++) + { + if (ffn_[nnIdx[doubleEdges[j].index]] != nnIdx[i]) + visibility = isVisible(uvn_nn[i], uvn_nn[doubleEdges[j].index], + doubleEdges[j].first, Eigen::Vector2f::Zero()); + if (!visibility) + break; + if (sfn_[nnIdx[doubleEdges[j].index]] != nnIdx[i]) + visibility = isVisible(uvn_nn[i], uvn_nn[doubleEdges[j].index], + doubleEdges[j].second, Eigen::Vector2f::Zero()); + if (!visibility == false) + break; + } + angles_[i].visible = visibility; + } + + // Selecting first two visible free neighbors + bool not_found = true; + PointId left = 1; + do + { + while ((left < nnn_) && + ((!angles_[left].visible) || stateSet(nnIdx[left]))) + left++; + if (left >= nnn_) + break; + else + { + int right = left+1; + do + { + while ((right < nnn_) && + ((!angles_[right].visible) || stateSet(nnIdx[right]))) + right++; + if (right >= nnn_) + break; + else if ((getCoord(nnIdx[left]) - + getCoord(nnIdx[right])).squaredNorm () > sqr_max_edge) + right++; + else + { + addFringePoint (nnIdx[right], R_); + addFringePoint (nnIdx[left], nnIdx[right]); + addFringePoint (R_, nnIdx[left]); + state_[R_] = state_[nnIdx[left]] = + state_[nnIdx[right]] = GP3Type::FRINGE; + ffn_[R_] = nnIdx[left]; + sfn_[R_] = nnIdx[right]; + ffn_[nnIdx[left]] = nnIdx[right]; + sfn_[nnIdx[left]] = R_; + ffn_[nnIdx[right]] = R_; + sfn_[nnIdx[right]] = nnIdx[left]; + addTriangle (R_, nnIdx[left], nnIdx[right]); + nr_parts++; + not_found = false; + break; + } + } + while (true); + left++; + } + } + while (not_found); + } + + // Set the index of the first free point in is_free. + done = true; + auto it = std::find(state_.begin(), state_.end(), GP3Type::FREE); + if (it == state_.end()) + done = true; + else + { + done = false; + is_free = std::distance(state_.begin(), it); + } + + is_fringe = true; + while (is_fringe) + { + is_fringe = false; + + int fqSize = static_cast (fringe_queue_.size ()); + while ((fqIdx < fqSize) && + (state_[fringe_queue_[fqIdx]] != GP3Type::FRINGE)) + fqIdx++; + + // an unfinished fringe point is found + if (fqIdx >= fqSize) + { + continue; + } + + R_ = fringe_queue_[fqIdx]; + is_fringe = true; + + if (ffn_[R_] == sfn_[R_]) + { + state_[R_] = GP3Type::COMPLETED; + continue; + } + tree.knnSearch(R_, nnn_, &nnIdx, &sqrDists); + +/** + // Search tree returns indices into the original cloud, but we are working with indices TODO: make that optional! + for (int i = 1; i < nnn_; i++) + { + //if (point2index[nnIdx[i]] == -1) + // std::cerr << R_ << " [" << indices_->at (R_) << "] " << i << ": " << nnIdx[i] << " / " << point2index[nnIdx[i]] << std::endl; + nnIdx[i] = point2index[nnIdx[i]]; + } +**/ + + Eigen::Vector3d coord(getCoord(R_)); + // Locating FFN and SFN to adapt distance threshold + double sqr_source_dist = (coord - getCoord[source_[R_]]).squaredNorm (); + double sqr_ffn_dist = (coord - getCoord[ffn_[R_]]).squaredNorm (); + double sqr_sfn_dist = (coord - getCoord[sfn_[R_]]).squaredNorm (); + double max_sqr_fn_dist = (std::max)(sqr_ffn_dist, sqr_sfn_dist); + double sqr_dist_threshold = (std::min)(sqr_max_edge, sqr_mu * sqrDists[1]); + if (max_sqr_fn_dist > sqrDists[nnn_-1]) + { + if (0 == increase_nnn4fn) + log()->get(LogLevel::Warning) << "Not enough neighbors are " + "considered: ffn or sfn out of range! Consider increasing " + "nnn_... Setting R = " << R_ << " to be BOUNDARY!\n"; + increase_nnn4fn++; + state_[R_] = GP3Type::BOUNDARY; + continue; + } + double max_sqr_fns_dist = (std::max)(sqr_source_dist, max_sqr_fn_dist); + if (max_sqr_fns_dist > sqrDists[nnn_-1]) + { + if (0 == increase_nnn4s) + log()->get(LogLevel::Warning) << "Not enough neighbors are " + "considered: source of R " << R_ << " is out of range! " + "Consider increasing nnn_...\n"; + increase_nnn4s++; + } + + // Get the normal estimate at the current point + const Eigen::Vector3d nc(getNormalCoord(R_)); + + // Get a coordinate system that lies on a plane defined by its normal + v_ = nc.unitOrthogonal (); + u_ = nc.cross (v_); + + const Eigen::Vector3d coord(getCoord(R_)); + // Projecting point onto the surface + float dist = nc.dot (coord); + proj_qp_ = coord - dist * nc; + + // Converting coords, calculating angles and saving the projected + // near boundary edges + int nr_edge = 0; + std::vector doubleEdges; + // nearest neighbor with index 0 is the query point R_ itself + //ABELL - Check this. + for (int i = 1; i < nnn_; i++) + { + tmp_ = coords_[nnIdx[i]] - proj_qp_; + uvn_nn[i][0] = tmp_.dot(u_); + uvn_nn[i][1] = tmp_.dot(v_); + + // Computing the angle between each neighboring point and the + // query point itself + angles_[i].angle = atan2(uvn_nn[i][1], uvn_nn[i][0]); + // initializing angle descriptors + angles_[i].index = nnIdx[i]; + angles_[i].nnIndex = i; + if ( + (state_[nnIdx[i]] == GP3Type::COMPLETED) || + (state_[nnIdx[i]] == GP3Type::BOUNDARY) || + (state_[nnIdx[i]] == GP3Type::NONE) || + (nnIdx[i] == -1) /// NOTE: discarding NaN points and those that are not in indices_ + || (sqrDists[i] > sqr_dist_threshold) + ) + angles_[i].visible = false; + else + angles_[i].visible = true; + if ((ffn_[R_] == nnIdx[i]) || (sfn_[R_] == nnIdx[i])) + angles_[i].visible = true; + bool same_side = true; + const Eigen::Vector3f neighbor_normal = getNormalCoord(nnIdx[i]); + double cosine = nc.dot (neighbor_normal); + if (cosine > 1) cosine = 1; + if (cosine < -1) cosine = -1; + double angle = acos (cosine); + if ((!consistent_) && (angle > M_PI/2)) + angle = M_PI - angle; + if (angle > eps_angle_) + { + angles_[i].visible = false; + same_side = false; + } + // Saving the edges between nearby boundary points + if ((i!=0) && (same_side) && + ((state_[nnIdx[i]] == GP3Type::FRINGE) || + (state_[nnIdx[i]] == GP3Type::BOUNDARY))) + { + doubleEdge e; + e.index = i; + nr_edge++; + tmp_ = getCoord[ffn_[nnIdx[i]]] - proj_qp_; + e.first[0] = tmp_.dot(u_); + e.first[1] = tmp_.dot(v_); + tmp_ = getCoord[sfn_[nnIdx[i]]] - proj_qp_; + e.second[0] = tmp_.dot(u_); + e.second[1] = tmp_.dot(v_); + doubleEdges.push_back(e); + // Pruning by visibility criterion + if ((state_[nnIdx[i]] == GP3Type::FRINGE) && + (ffn_[R_] != nnIdx[i]) && (sfn_[R_] != nnIdx[i])) + { + double angle1 = atan2(e.first[1] - uvn_nn[i][1], + e.first[0] - uvn_nn[i][0]); + double angle2 = atan2(e.second[1] - uvn_nn[i][1], + e.second[0] - uvn_nn[i][0]); + double angleMin, angleMax; + if (angle1 < angle2) + { + angleMin = angle1; + angleMax = angle2; + } + else + { + angleMin = angle2; + angleMax = angle1; + } + double angleR = angles_[i].angle + M_PI; + if (angleR >= 2*M_PI) + angleR -= 2*M_PI; + if ((source_[nnIdx[i]] == ffn_[nnIdx[i]]) || (source_[nnIdx[i]] == sfn_[nnIdx[i]])) + { + if ((angleMax - angleMin) < M_PI) + { + if ((angleMin < angleR) && (angleR < angleMax)) + angles_[i].visible = false; + } + else + { + if ((angleR < angleMin) || (angleMax < angleR)) + angles_[i].visible = false; + } + } + else + { + tmp_ = getCoord(source_[nnIdx[i]]) - proj_qp_; + uvn_s[0] = tmp_.dot(u_); + uvn_s[1] = tmp_.dot(v_); + double angleS = atan2(uvn_s[1] - uvn_nn[i][1], + uvn_s[0] - uvn_nn[i][0]); + if ((angleMin < angleS) && (angleS < angleMax)) + { + if ((angleMin < angleR) && (angleR < angleMax)) + angles_[i].visible = false; + } + else + { + if ((angleR < angleMin) || (angleMax < angleR)) + angles_[i].visible = false; + } + } + } + } + } + angles_[0].visible = false; + + // Verify the visibility of each potential new vertex + // nearest neighbor with index 0 is the query point R_ itself + //ABELL - Check above (AGAIN?!) + for (int i = 1; i < nnn_; i++) + if ((angles_[i].visible) && + (ffn_[R_] != nnIdx[i]) && (sfn_[R_] != nnIdx[i])) + { + bool visibility = true; + for (int j = 0; j < nr_edge; j++) + { + if (doubleEdges[j].index != i) + { + int f = ffn_[nnIdx[doubleEdges[j].index]]; + if ((f != nnIdx[i]) && (f != R_)) + visibility = isVisible(uvn_nn[i], uvn_nn[doubleEdges[j].index], + doubleEdges[j].first, Eigen::Vector2f::Zero()); + if (visibility == false) + break; + + int s = sfn_[nnIdx[doubleEdges[j].index]]; + if ((s != nnIdx[i]) && (s != R_)) + visibility = isVisible(uvn_nn[i], uvn_nn[doubleEdges[j].index], + doubleEdges[j].second, Eigen::Vector2f::Zero()); + if (visibility == false) + break; + } + } + angles_[i].visible = visibility; + } + + // Sorting angles + auto comp = [](const nnAngle& a1, const nnAngle + std::sort(angles_.begin (), angles_.end (), + [](const nnAngle& a1, const nnAngle& a2) + { + return (a1.visible == a2.visible ? + a1.angle < a2.angle : + a1.visible); + }); + + // Triangulating + if (angles_[2].visible == false) + { + if ( !( (angles_[0].index == ffn_[R_] && angles_[1].index == sfn_[R_]) + || (angles_[0].index == sfn_[R_] && angles_[1].index == ffn_[R_]) )) + { + state_[R_] = GP3Type::BOUNDARY; + } + else + { + if ((source_[R_] == angles_[0].index) || + (source_[R_] == angles_[1].index)) + state_[R_] = GP3Type::BOUNDARY; + else + { + if (sqr_max_edge < + (coords_[ffn_[R_]] - coords_[sfn_[R_]]).squaredNorm ()) + { + state_[R_] = GP3Type::BOUNDARY; + } + else + { + tmp_ = coords_[source_[R_]] - proj_qp_; + uvn_s[0] = tmp_.dot(u_); + uvn_s[1] = tmp_.dot(v_); + double angleS = atan2(uvn_s[1], uvn_s[0]); + double dif = angles_[1].angle - angles_[0].angle; + if ((angles_[0].angle < angleS) && (angleS < angles_[1].angle)) + { + if (dif < 2*M_PI - maximum_angle_) + state_[R_] = GP3Type::BOUNDARY; + else + closeTriangle(); + } + else + { + if (dif >= maximum_angle_) + state_[R_] = GP3Type::BOUNDARY; + else + closeTriangle(); + } + } + } + } + continue; + } + + // Finding the FFN and SFN + int start = -1, end = -1; + for (int i=0; i end)) + need_invert = true; + } + else + { + if (nCB != -1) + { + if ((nCB == start) || (nCB == end)) + { + bool inside_CB = false; + bool outside_CB = false; + for (int i=0; i angles_[start].angle) && + (angles_[nCB].angle < angles_[end].angle)) + need_invert = true; + } + } + else + { + if (start == end-1) + need_invert = true; + } + } + } + else if ((angles_[start].angle < angles_[sourceIdx].angle) && + (angles_[sourceIdx].angle < angles_[end].angle)) + need_invert = true; + } + + // switching start and end if necessary + if (need_invert) + std::swap(start, end); + + // Arranging visible nnAngles in the order they need to be connected and + // compute the maximal angle difference between two consecutive visible + // angles + bool is_boundary = false, is_skinny = false; + std::vector gaps (nnn_, false); + std::vector skinny (nnn_, false); + std::vector dif (nnn_); + std::vector angleIdx; angleIdx.reserve (nnn_); + if (start > end) + { + for (int j=start; j::iterator first_gap_after = angleIdx.end (); + std::vector::iterator last_gap_before = angleIdx.begin (); + int nr_gaps = 0; + for (auto it = angleIdx.begin (); it != angleIdx.end () - 1; it++) + { + if (gaps[*it]) + { + nr_gaps++; + if (first_gap_after == angleIdx.end()) + first_gap_after = it; + last_gap_before = it + 1; + } + } + if (nr_gaps > 1) + { + angleIdx.erase(first_gap_after+1, last_gap_before); + } + + // Neglecting points that would form skinny triangles (if possible) + if (is_skinny) + { + double angle_so_far = 0, angle_would_be; + double max_combined_angle = + (std::min)(maximum_angle_, M_PI-2*minimum_angle_); + Eigen::Vector2f X; + Eigen::Vector2f S1; + Eigen::Vector2f S2; + std::vector to_erase; + for (auto it = angleIdx.begin()+1; it != angleIdx.end()-1; it++) + { + if (gaps[*(it-1)]) + angle_so_far = 0; + else + angle_so_far += dif[*(it-1)]; + if (gaps[*it]) + angle_would_be = angle_so_far; + else + angle_would_be = angle_so_far + dif[*it]; + if ( + (skinny[*it] || skinny[*(it-1)]) && + (!stateSet(angles_[*it].index) || + !stateSet(angles_[*(it-1)].index)) && + ((!gaps[*it]) || + (angles_[*it].nnIndex > angles_[*(it-1)].nnIndex)) && + ((!gaps[*(it-1)]) || + (angles_[*it].nnIndex > angles_[*(it+1)].nnIndex)) && + (angle_would_be < max_combined_angle) + ) + { + if (gaps[*(it-1)]) + { + gaps[*it] = true; + to_erase.push_back(*it); + } + else if (gaps[*it]) + { + gaps[*(it-1)] = true; + to_erase.push_back(*it); + } + else + { + int erased_idx = static_cast (to_erase.size ()) -1; + for (auto prev_it = it-1; + (erased_idx != -1) && (it != angleIdx.begin()); it--) + if (*it == to_erase[erased_idx]) + erased_idx--; + else + break; + bool can_delete = true; + for (auto curr_it = prev_it+1; curr_it != it+1; curr_it++) + { + tmp_ = getCoord(angles_[*curr_it].index) - proj_qp_; + X[0] = tmp_.dot(u_); + X[1] = tmp_.dot(v_); + tmp_ = getCoord(angles_[*prev_it].index) - proj_qp_; + S1[0] = tmp_.dot(u_); + S1[1] = tmp_.dot(v_); + tmp_ = getCoord(angles_[*(it+1)].index) - proj_qp_; + S2[0] = tmp_.dot(u_); + S2[1] = tmp_.dot(v_); + // check for inclusions + if (isVisible(X,S1,S2)) + { + can_delete = false; + angle_so_far = 0; + break; + } + } + if (can_delete) + { + to_erase.push_back(*it); + } + } + } + else + angle_so_far = 0; + } + for (auto it = to_erase.begin(); it != to_erase.end(); it++) + { + for (auto iter = angleIdx.begin(); iter != angleIdx.end(); iter++) + if (*it == *iter) + { + angleIdx.erase(iter); + break; + } + } + } + + // Writing edges and updating edge-front + changed_1st_fn_ = false; + changed_2nd_fn_ = false; +//ABELL - This seems bad - sentinel for index that's positive. + new2boundary_ = -1; + for (auto it = angleIdx.begin() + 1; it != angleIdx.end() - 1; it++) + { + current_index_ = angles_[*it].index; + + is_current_free_ = false; + if (!stateSet(current_index_)) + { + state_[current_index_] = GP3Type::FRINGE; + is_current_free_ = true; + } + else if (!already_connected_) + { + prev_is_ffn_ = (ffn_[current_index_] == + angles_[*(it-1)].index) && (!gaps[*(it-1)]); + prev_is_sfn_ = (sfn_[current_index_] == + angles_[*(it-1)].index) && (!gaps[*(it-1)]); + next_is_ffn_ = (ffn_[current_index_] == + angles_[*(it+1)].index) && (!gaps[*it]); + next_is_sfn_ = (sfn_[current_index_] == + angles_[*(it+1)].index) && (!gaps[*it]); + if (!prev_is_ffn_ && !next_is_sfn_ && !prev_is_sfn_ && !next_is_ffn_) + { + nr_touched++; + } + } + + if (gaps[*it]) + if (gaps[*(it-1)]) + { + if (is_current_free_) + state_[current_index_] = GP3Type::NONE; /// TODO: document! + } + + else // (gaps[*it]) && ^(gaps[*(it-1)]) + { + addTriangle (current_index_, angles_[*(it-1)].index, R_); + addFringePoint (current_index_, R_); + new2boundary_ = current_index_; + if (!already_connected_) + connectPoint (angles_[*(it-1)].index, R_, angles_[*(it+1)].index, + uvn_nn[angles_[*it].nnIndex], + uvn_nn[angles_[*(it-1)].nnIndex], uvn_nn_qp_zero); + else + already_connected_ = false; + if (ffn_[R_] == angles_[*(angleIdx.begin())].index) + { + ffn_[R_] = new2boundary_; + } + else if (sfn_[R_] == angles_[*(angleIdx.begin())].index) + { + sfn_[R_] = new2boundary_; + } + } + + else // ^(gaps[*it]) + if (gaps[*(it-1)]) + { + addFringePoint (current_index_, R_); + new2boundary_ = current_index_; + if (!already_connected_) + connectPoint(R_, angles_[*(it+1)].index, + (it+2) == angleIdx.end() ? -1 : angles_[*(it+2)].index, + uvn_nn[angles_[*it].nnIndex], uvn_nn_qp_zero, + uvn_nn[angles_[*(it+1)].nnIndex]); + else + already_connected_ = false; + if (ffn_[R_] == angles_[*(angleIdx.end()-1)].index) + { + ffn_[R_] = new2boundary_; + } + else if (sfn_[R_] == angles_[*(angleIdx.end()-1)].index) + { + sfn_[R_] = new2boundary_; + } + } + + else // ^(gaps[*it]) && ^(gaps[*(it-1)]) + { + addTriangle (current_index_, angles_[*(it-1)].index, R_); + addFringePoint (current_index_, R_); + if (!already_connected_) + connectPoint (angles_[*(it-1)].index, angles_[*(it+1)].index, + (it+2) == angleIdx.end() ? -1 : + gaps[*(it+1)] ? R_ : + angles_[*(it+2)].index, + uvn_nn[angles_[*it].nnIndex], + uvn_nn[angles_[*(it-1)].nnIndex], + uvn_nn[angles_[*(it+1)].nnIndex]); + else + already_connected_ = false; + } + } + + // Finishing up R_ + if (ffn_[R_] == sfn_[R_]) + { + state_[R_] = GP3Type::COMPLETED; + } + if (!gaps[*(angleIdx.end()-2)]) + { + addTriangle (angles_[*(angleIdx.end()-2)].index, + angles_[*(angleIdx.end()-1)].index, R_); + addFringePoint (angles_[*(angleIdx.end()-2)].index, R_); + if (R_ == ffn_[angles_[*(angleIdx.end()-1)].index]) + { + if (angles_[*(angleIdx.end()-2)].index == + sfn_[angles_[*(angleIdx.end()-1)].index]) + { + state_[angles_[*(angleIdx.end()-1)].index] = GP3Type::COMPLETED; + } + else + { + ffn_[angles_[*(angleIdx.end()-1)].index] = + angles_[*(angleIdx.end()-2)].index; + } + } + else if (R_ == sfn_[angles_[*(angleIdx.end()-1)].index]) + { + if (angles_[*(angleIdx.end()-2)].index == + ffn_[angles_[*(angleIdx.end()-1)].index]) + { + state_[angles_[*(angleIdx.end()-1)].index] = GP3Type::COMPLETED; + } + else + { + sfn_[angles_[*(angleIdx.end()-1)].index] = + angles_[*(angleIdx.end()-2)].index; + } + } + } + if (!gaps[*(angleIdx.begin())]) + { + if (R_ == ffn_[angles_[*(angleIdx.begin())].index]) + { + if (angles_[*(angleIdx.begin()+1)].index == + sfn_[angles_[*(angleIdx.begin())].index]) + { + state_[angles_[*(angleIdx.begin())].index] = GP3Type::COMPLETED; + } + else + { + ffn_[angles_[*(angleIdx.begin())].index] = + angles_[*(angleIdx.begin()+1)].index; + } + } + else if (R_ == sfn_[angles_[*(angleIdx.begin())].index]) + { + if (angles_[*(angleIdx.begin()+1)].index == + ffn_[angles_[*(angleIdx.begin())].index]) + { + state_[angles_[*(angleIdx.begin())].index] = GP3Type::COMPLETED; + } + else + { + sfn_[angles_[*(angleIdx.begin())].index] = + angles_[*(angleIdx.begin()+1)].index; + } + } + } + } + } + log()->get(LogLevel::Debug) << "Number of triangles: " << grid_.size() << + ".\n"); + log()->get(LogLevel::Debug) << "Number of unconnected parts: " << nr_parts << + ".\n"); + if (increase_nnn4fn > 0) + log()->get(LogLevel::Warning) << "Number of neighborhood size " + "increase requests for fringe neighbors: " << increase_nnn4fn << + ".\n"); + if (increase_nnn4s > 0) + log()->get(LogLevel::Warning) << "Number of neighborhood size " + "increase requests for source: " << increase_nnn4s << ".\n"; + if (increase_dist > 0) + log()->get(LogLevel::Warning) << "Number of automatic maximum " + "distance increases: " << increase_dist << ".\n"; + + // sorting and removing doubles from fringe queue + std::sort (fringe_queue_.begin(), fringe_queue_.end ()); + fringe_queue_.erase (std::unique(fringe_queue_.begin(), fringe_queue_.end()), + fringe_queue_.end ()); + log()->get(LogLevel::Debug) << "Number of processed points: " << + fringe_queue_.size() << " / " << view.size() << "!\n"; + view_ = nullptr; + return (true); +} + +GreedyProjection::closeTriangle () +{ + state_[R_] = GP3Type::COMPLETED; + addTriangle (angles_[0].index, angles_[1].index, R_); + for (int aIdx=0; aIdx<2; aIdx++) + { + if (ffn_[angles_[aIdx].index] == R_) + { + if (sfn_[angles_[aIdx].index] == angles_[(aIdx+1)%2].index) + { + state_[angles_[aIdx].index] = GP3Type::COMPLETED; + } + else + { + ffn_[angles_[aIdx].index] = angles_[(aIdx+1)%2].index; + } + } + else if (sfn_[angles_[aIdx].index] == R_) + { + if (ffn_[angles_[aIdx].index] == angles_[(aIdx+1)%2].index) + { + state_[angles_[aIdx].index] = GP3Type::COMPLETED; + } + else + { + sfn_[angles_[aIdx].index] = angles_[(aIdx+1)%2].index; + } + } + } +} + +void GreedyProjection::connectPoint ( + const int prev_index, const int next_index, const int next_next_index, + const Eigen::Vector2d &uvn_current, + const Eigen::Vector2d &uvn_prev, + const Eigen::Vector2d &uvn_next) +{ + if (is_current_free_) + { + ffn_[current_index_] = prev_index; + sfn_[current_index_] = next_index; + } + else + { + if ((prev_is_ffn_ && next_is_sfn_) || (prev_is_sfn_ && next_is_ffn_)) + state_[current_index_] = GP3Type::COMPLETED; + else if (prev_is_ffn_ && !next_is_sfn_) + ffn_[current_index_] = next_index; + else if (next_is_ffn_ && !prev_is_sfn_) + ffn_[current_index_] = prev_index; + else if (prev_is_sfn_ && !next_is_ffn_) + sfn_[current_index_] = next_index; + else if (next_is_sfn_ && !prev_is_ffn_) + sfn_[current_index_] = prev_index; + else + { + bool found_triangle = false; + if ((prev_index != R_) && + ((ffn_[current_index_] == ffn_[prev_index]) || + (ffn_[current_index_] == sfn_[prev_index]))) + { + found_triangle = true; + addTriangle (current_index_, ffn_[current_index_], prev_index); + state_[prev_index] = GP3Type::COMPLETED; + state_[ffn_[current_index_]] = GP3Type::COMPLETED; + ffn_[current_index_] = next_index; + } + else if ((prev_index != R_) && + ((sfn_[current_index_] == ffn_[prev_index]) || + (sfn_[current_index_] == sfn_[prev_index]))) + { + found_triangle = true; + addTriangle (current_index_, sfn_[current_index_], prev_index); + state_[prev_index] = GP3Type::COMPLETED; + state_[sfn_[current_index_]] = GP3Type::COMPLETED; + sfn_[current_index_] = next_index; + } + else if (stateSet(next_index)) + { + if ((ffn_[current_index_] == ffn_[next_index]) || + (ffn_[current_index_] == sfn_[next_index])) + { + found_triangle = true; + addTriangle (current_index_, ffn_[current_index_], next_index); + + if (ffn_[current_index_] == ffn_[next_index]) + { + ffn_[next_index] = current_index_; + } + else + { + sfn_[next_index] = current_index_; + } + state_[ffn_[current_index_]] = GP3Type::COMPLETED; + ffn_[current_index_] = prev_index; + } + else if ((sfn_[current_index_] == ffn_[next_index]) || + (sfn_[current_index_] == sfn_[next_index])) + { + found_triangle = true; + addTriangle (current_index_, sfn_[current_index_], next_index); + + if (sfn_[current_index_] == ffn_[next_index]) + { + ffn_[next_index] = current_index_; + } + else + { + sfn_[next_index] = current_index_; + } + state_[sfn_[current_index_]] = GP3Type::COMPLETED; + sfn_[current_index_] = prev_index; + } + } + + if (found_triangle) + { + } + else + { + tmp_ = getCoord(ffn_[current_index_]) - proj_qp_; + uvn_ffn_[0] = tmp_.dot(u_); + uvn_ffn_[1] = tmp_.dot(v_); + tmp_ = getCoord(sfn_[current_index_]) - proj_qp_; + uvn_sfn_[0] = tmp_.dot(u_); + uvn_sfn_[1] = tmp_.dot(v_); + bool prev_ffn = isVisible(uvn_prev, uvn_next, uvn_current, uvn_ffn_) && + isVisible(uvn_prev, uvn_sfn_, uvn_current, uvn_ffn_); + bool prev_sfn = isVisible(uvn_prev, uvn_next, uvn_current, uvn_sfn_) && + isVisible(uvn_prev, uvn_ffn_, uvn_current, uvn_sfn_); + bool next_ffn = isVisible(uvn_next, uvn_prev, uvn_current, uvn_ffn_) && + isVisible(uvn_next, uvn_sfn_, uvn_current, uvn_ffn_); + bool next_sfn = isVisible(uvn_next, uvn_prev, uvn_current, uvn_sfn_) && + isVisible(uvn_next, uvn_ffn_, uvn_current, uvn_sfn_); + int min_dist = -1; + if (prev_ffn && next_sfn && prev_sfn && next_ffn) + { + /* should be never the case */ + double prev2f = (getCoord[ffn_[current_index_]] - + getCoord[prev_index]).squaredNorm (); + double next2s = (getCoord[sfn_[current_index_]] - + getCoord[next_index]).squaredNorm (); + double prev2s = (getCoord[sfn_[current_index_]] - + getCoord[prev_index]).squaredNorm (); + double next2f = (getCoord[ffn_[current_index_]] - + getCoord[next_index]).squaredNorm (); + if (prev2f < prev2s) + { + if (prev2f < next2f) + { + if (prev2f < next2s) + min_dist = 0; + else + min_dist = 3; + } + else + { + if (next2f < next2s) + min_dist = 2; + else + min_dist = 3; + } + } + else + { + if (prev2s < next2f) + { + if (prev2s < next2s) + min_dist = 1; + else + min_dist = 3; + } + else + { + if (next2f < next2s) + min_dist = 2; + else + min_dist = 3; + } + } + } + else if (prev_ffn && next_sfn) + { + /* a clear case */ + double prev2f = (getCoord[ffn_[current_index_]] - + getCoord[prev_index]).squaredNorm (); + double next2s = (getCoord[sfn_[current_index_]] - + getCoord[next_index]).squaredNorm (); + if (prev2f < next2s) + min_dist = 0; + else + min_dist = 3; + } + else if (prev_sfn && next_ffn) + { + /* a clear case */ + double prev2s = (getCoord[sfn_[current_index_]] - + getCoord[prev_index]).squaredNorm (); + double next2f = (getCoord[ffn_[current_index_]] - + getCoord[next_index]).squaredNorm (); + if (prev2s < next2f) + min_dist = 1; + else + min_dist = 2; + } + /* straightforward cases */ + else if (prev_ffn && !next_sfn && !prev_sfn && !next_ffn) + min_dist = 0; + else if (!prev_ffn && !next_sfn && prev_sfn && !next_ffn) + min_dist = 1; + else if (!prev_ffn && !next_sfn && !prev_sfn && next_ffn) + min_dist = 2; + else if (!prev_ffn && next_sfn && !prev_sfn && !next_ffn) + min_dist = 3; + /* messed up cases */ + else if (prev_ffn) + { + double prev2f = (getCoord[ffn_[current_index_]] - + getCoord[prev_index]).squaredNorm (); + if (prev_sfn) + { + double prev2s = (getCoord[sfn_[current_index_]] - + getCoord[prev_index]).squaredNorm (); + if (prev2s < prev2f) + min_dist = 1; + else + min_dist = 0; + } + else if (next_ffn) + { + double next2f = (getCoord[ffn_[current_index_]] - + getCoord[next_index]).squaredNorm (); + if (next2f < prev2f) + min_dist = 2; + else + min_dist = 0; + } + } + else if (next_sfn) + { + double next2s = (getCoord[sfn_[current_index_]] - + getCoord[next_index]).squaredNorm (); + if (prev_sfn) + { + double prev2s = (getCoord[sfn_[current_index_]] - + getCoord[prev_index]).squaredNorm (); + if (prev2s < next2s) + min_dist = 1; + else + min_dist = 3; + } + else if (next_ffn) + { + double next2f = (getCoord[ffn_[current_index_]] - + getCoord[next_index]).squaredNorm (); + if (next2f < next2s) + min_dist = 2; + else + min_dist = 3; + } + } + switch (min_dist) + { + case 0://prev2f: + { + addTriangle (current_index_, ffn_[current_index_], prev_index); + + /* updating prev_index */ + if (ffn_[prev_index] == current_index_) + { + ffn_[prev_index] = ffn_[current_index_]; + } + else if (sfn_[prev_index] == current_index_) + { + sfn_[prev_index] = ffn_[current_index_]; + } + else if (ffn_[prev_index] == R_) + { + changed_1st_fn_ = true; + ffn_[prev_index] = ffn_[current_index_]; + } + else if (sfn_[prev_index] == R_) + { + changed_1st_fn_ = true; + sfn_[prev_index] = ffn_[current_index_]; + } + else if (prev_index == R_) + { + new2boundary_ = ffn_[current_index_]; + } + + /* updating ffn */ + if (ffn_[ffn_[current_index_]] == current_index_) + { + ffn_[ffn_[current_index_]] = prev_index; + } + else if (sfn_[ffn_[current_index_]] == current_index_) + { + sfn_[ffn_[current_index_]] = prev_index; + } + + /* updating current */ + ffn_[current_index_] = next_index; + + break; + } + case 1://prev2s: + { + addTriangle (current_index_, sfn_[current_index_], prev_index); + + /* updating prev_index */ + if (ffn_[prev_index] == current_index_) + { + ffn_[prev_index] = sfn_[current_index_]; + } + else if (sfn_[prev_index] == current_index_) + { + sfn_[prev_index] = sfn_[current_index_]; + } + else if (ffn_[prev_index] == R_) + { + changed_1st_fn_ = true; + ffn_[prev_index] = sfn_[current_index_]; + } + else if (sfn_[prev_index] == R_) + { + changed_1st_fn_ = true; + sfn_[prev_index] = sfn_[current_index_]; + } + else if (prev_index == R_) + { + new2boundary_ = sfn_[current_index_]; + } + + /* updating sfn */ + if (ffn_[sfn_[current_index_]] == current_index_) + { + ffn_[sfn_[current_index_]] = prev_index; + } + else if (sfn_[sfn_[current_index_]] == current_index_) + { + sfn_[sfn_[current_index_]] = prev_index; + } + + /* updating current */ + sfn_[current_index_] = next_index; + + break; + } + case 2://next2f: + { + addTriangle (current_index_, ffn_[current_index_], next_index); + int neighbor_update = next_index; + + /* updating next_index */ + if (!stateSet(next_index)) + { + state_[next_index] = GP3Type::FRINGE; + ffn_[next_index] = current_index_; + sfn_[next_index] = ffn_[current_index_]; + } + else + { + if (ffn_[next_index] == R_) + { + changed_2nd_fn_ = true; + ffn_[next_index] = ffn_[current_index_]; + } + else if (sfn_[next_index] == R_) + { + changed_2nd_fn_ = true; + sfn_[next_index] = ffn_[current_index_]; + } + else if (next_index == R_) + { + new2boundary_ = ffn_[current_index_]; + if (next_next_index == new2boundary_) + already_connected_ = true; + } + else if (ffn_[next_index] == next_next_index) + { + already_connected_ = true; + ffn_[next_index] = ffn_[current_index_]; + } + else if (sfn_[next_index] == next_next_index) + { + already_connected_ = true; + sfn_[next_index] = ffn_[current_index_]; + } + else + { + tmp_ = coords_[ffn_[next_index]] - proj_qp_; + uvn_next_ffn_[0] = tmp_.dot(u_); + uvn_next_ffn_[1] = tmp_.dot(v_); + tmp_ = coords_[sfn_[next_index]] - proj_qp_; + uvn_next_sfn_[0] = tmp_.dot(u_); + uvn_next_sfn_[1] = tmp_.dot(v_); + + bool ffn_next_ffn = + isVisible(uvn_next_ffn_, uvn_next, uvn_current, uvn_ffn_) && + isVisible(uvn_next_ffn_, uvn_next, uvn_next_sfn_, uvn_ffn_); + bool sfn_next_ffn = + isVisible(uvn_next_sfn_, uvn_next, uvn_current, uvn_ffn_) && + isVisible(uvn_next_sfn_, uvn_next, uvn_next_ffn_, uvn_ffn_); + + int connect2ffn = -1; + if (ffn_next_ffn && sfn_next_ffn) + { + double fn2f = (getCoord[ffn_[current_index_]] - + getCoord[ffn_[next_index]]).squaredNorm (); + double sn2f = (getCoord[ffn_[current_index_]] - + getCoord[sfn_[next_index]]).squaredNorm (); + if (fn2f < sn2f) + connect2ffn = 0; + else + connect2ffn = 1; + } + else if (ffn_next_ffn) + connect2ffn = 0; + else if (sfn_next_ffn) + connect2ffn = 1; + + switch (connect2ffn) + { + case 0: // ffn[next] + { + addTriangle (next_index, ffn_[current_index_], + ffn_[next_index]); + neighbor_update = ffn_[next_index]; + + /* ffn[next_index] */ + if ((ffn_[ffn_[next_index]] == ffn_[current_index_]) || + (sfn_[ffn_[next_index]] == ffn_[current_index_])) + { + state_[ffn_[next_index]] = GP3Type::COMPLETED; + } + else if (ffn_[ffn_[next_index]] == next_index) + { + ffn_[ffn_[next_index]] = ffn_[current_index_]; + } + else if (sfn_[ffn_[next_index]] == next_index) + { + sfn_[ffn_[next_index]] = ffn_[current_index_]; + } + + ffn_[next_index] = current_index_; + + break; + } + case 1: // sfn[next] + { + addTriangle (next_index, ffn_[current_index_], + sfn_[next_index]); + neighbor_update = sfn_[next_index]; + + /* sfn[next_index] */ + if ((ffn_[sfn_[next_index]] = ffn_[current_index_]) || + (sfn_[sfn_[next_index]] == ffn_[current_index_])) + { + state_[sfn_[next_index]] = GP3Type::COMPLETED; + } + else if (ffn_[sfn_[next_index]] == next_index) + { + ffn_[sfn_[next_index]] = ffn_[current_index_]; + } + else if (sfn_[sfn_[next_index]] == next_index) + { + sfn_[sfn_[next_index]] = ffn_[current_index_]; + } + + sfn_[next_index] = current_index_; + + break; + } + default:; + } + } + } + + /* updating ffn */ + if ((ffn_[ffn_[current_index_]] == neighbor_update) || + (sfn_[ffn_[current_index_]] == neighbor_update)) + { + state_[ffn_[current_index_]] = GP3Type::COMPLETED; + } + else if (ffn_[ffn_[current_index_]] == current_index_) + { + ffn_[ffn_[current_index_]] = neighbor_update; + } + else if (sfn_[ffn_[current_index_]] == current_index_) + { + sfn_[ffn_[current_index_]] = neighbor_update; + } + + /* updating current */ + ffn_[current_index_] = prev_index; + + break; + } + case 3://next2s: + { + addTriangle (current_index_, sfn_[current_index_], next_index); + int neighbor_update = next_index; + + /* updating next_index */ + if (!stateSet(next_index)) + { + state_[next_index] = GP3Type::FRINGE; + ffn_[next_index] = current_index_; + sfn_[next_index] = sfn_[current_index_]; + } + else + { + if (ffn_[next_index] == R_) + { + changed_2nd_fn_ = true; + ffn_[next_index] = sfn_[current_index_]; + } + else if (sfn_[next_index] == R_) + { + changed_2nd_fn_ = true; + sfn_[next_index] = sfn_[current_index_]; + } + else if (next_index == R_) + { + new2boundary_ = sfn_[current_index_]; + if (next_next_index == new2boundary_) + already_connected_ = true; + } + else if (ffn_[next_index] == next_next_index) + { + already_connected_ = true; + ffn_[next_index] = sfn_[current_index_]; + } + else if (sfn_[next_index] == next_next_index) + { + already_connected_ = true; + sfn_[next_index] = sfn_[current_index_]; + } + else + { + tmp_ = coords_[ffn_[next_index]] - proj_qp_; + uvn_next_ffn_[0] = tmp_.dot(u_); + uvn_next_ffn_[1] = tmp_.dot(v_); + tmp_ = coords_[sfn_[next_index]] - proj_qp_; + uvn_next_sfn_[0] = tmp_.dot(u_); + uvn_next_sfn_[1] = tmp_.dot(v_); + + bool ffn_next_sfn = + isVisible(uvn_next_ffn_, uvn_next, uvn_current, uvn_sfn_) && + isVisible(uvn_next_ffn_, uvn_next, uvn_next_sfn_, uvn_sfn_); + bool sfn_next_sfn = + isVisible(uvn_next_sfn_, uvn_next, uvn_current, uvn_sfn_) && + isVisible(uvn_next_sfn_, uvn_next, uvn_next_ffn_, uvn_sfn_); + + int connect2sfn = -1; + if (ffn_next_sfn && sfn_next_sfn) + { + double fn2s = (getCoord[sfn_[current_index_]] - + getCoord[ffn_[next_index]]).squaredNorm (); + double sn2s = (getCoord[sfn_[current_index_]] - + getCoord[sfn_[next_index]]).squaredNorm (); + if (fn2s < sn2s) + connect2sfn = 0; + else + connect2sfn = 1; + } + else if (ffn_next_sfn) + connect2sfn = 0; + else if (sfn_next_sfn) + connect2sfn = 1; + + switch (connect2sfn) + { + case 0: // ffn[next] + { + addTriangle (next_index, sfn_[current_index_], + ffn_[next_index]); + neighbor_update = ffn_[next_index]; + + /* ffn[next_index] */ + if ((ffn_[ffn_[next_index]] == sfn_[current_index_]) || + (sfn_[ffn_[next_index]] == sfn_[current_index_])) + { + state_[ffn_[next_index]] = GP3Type::COMPLETED; + } + else if (ffn_[ffn_[next_index]] == next_index) + { + ffn_[ffn_[next_index]] = sfn_[current_index_]; + } + else if (sfn_[ffn_[next_index]] == next_index) + { + sfn_[ffn_[next_index]] = sfn_[current_index_]; + } + + ffn_[next_index] = current_index_; + + break; + } + case 1: // sfn[next] + { + addTriangle (next_index, sfn_[current_index_], + sfn_[next_index]); + neighbor_update = sfn_[next_index]; + + /* sfn[next_index] */ + if ((ffn_[sfn_[next_index]] == sfn_[current_index_]) || + (sfn_[sfn_[next_index]] == sfn_[current_index_])) + { + state_[sfn_[next_index]] = GP3Type::COMPLETED; + } + else if (ffn_[sfn_[next_index]] == next_index) + { + ffn_[sfn_[next_index]] = sfn_[current_index_]; + } + else if (sfn_[sfn_[next_index]] == next_index) + { + sfn_[sfn_[next_index]] = sfn_[current_index_]; + } + + sfn_[next_index] = current_index_; + + break; + } + default:; + } + } + } + + /* updating sfn */ + if ((ffn_[sfn_[current_index_]] == neighbor_update) || + (sfn_[sfn_[current_index_]] == neighbor_update)) + { + state_[sfn_[current_index_]] = GP3Type::COMPLETED; + } + else if (ffn_[sfn_[current_index_]] == current_index_) + { + ffn_[sfn_[current_index_]] = neighbor_update; + } + else if (sfn_[sfn_[current_index_]] == current_index_) + { + sfn_[sfn_[current_index_]] = neighbor_update; + } + + sfn_[current_index_] = prev_index; + + break; + } + default:; + } + } + } + } +} + +/** +std::vector > +GreedyProjection::getTriangleList (const pcl::PolygonMesh &input) +{ + std::vector > triangleList (input.cloud.width * input.cloud.height); + + for (size_t i=0; i < input.polygons.size (); ++i) + for (size_t j=0; j < input.polygons[i].vertices.size (); ++j) + triangleList[input.polygons[i].vertices[j]].push_back (i); + return (triangleList); +} +**/ + +} // namespace pdal diff --git a/filters/GreedyProjection.hpp b/filters/GreedyProjection.hpp new file mode 100644 index 0000000000..79dc8a2fad --- /dev/null +++ b/filters/GreedyProjection.hpp @@ -0,0 +1,459 @@ +/* + * Software License Agreement (BSD License) + * + * Point Cloud Library (PCL) - www.pointclouds.org + * Copyright (c) 2010-2011, Willow Garage, Inc. + * + * 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 Willow Garage, 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. + * + * $Id$ + * + */ + +#pragma once + +#include +#include +#include +#include + +namespace pdal +{ + /** \brief Returns if a point X is visible from point R (or the origin) + * when taking into account the segment between the points S1 and S2 + * \param X 2D coordinate of the point + * \param S1 2D coordinate of the segment's first point + * \param S2 2D coordinate of the segment's secont point + * \param R 2D coorddinate of the reference point (defaults to 0,0) + * \ingroup surface + */ + inline bool + isVisible (const Eigen::Vector2d &X, const Eigen::Vector2d &S1, + const Eigen::Vector2d &S2, + const Eigen::Vector2d &R = Eigen::Vector2f::Zero ()) + { + double a0 = S1[1] - S2[1]; + double b0 = S2[0] - S1[0]; + double c0 = S1[0]*S2[1] - S2[0]*S1[1]; + double a1 = -X[1]; + double b1 = X[0]; + double c1 = 0; + if (R != Eigen::Vector2f::Zero()) + { + a1 += R[1]; + b1 -= R[0]; + c1 = R[0]*X[1] - X[0]*R[1]; + } + double div = a0*b1 - b0*a1; + double x = (b0*c1 - b1*c0) / div; + double y = (a1*c0 - a0*c1) / div; + + bool intersection_outside_XR; + if (R == Eigen::Vector2f::Zero()) + { + if (X[0] > 0) + intersection_outside_XR = (x <= 0) || (x >= X[0]); + else if (X[0] < 0) + intersection_outside_XR = (x >= 0) || (x <= X[0]); + else if (X[1] > 0) + intersection_outside_XR = (y <= 0) || (y >= X[1]); + else if (X[1] < 0) + intersection_outside_XR = (y >= 0) || (y <= X[1]); + else + intersection_outside_XR = true; + } + else + { + if (X[0] > R[0]) + intersection_outside_XR = (x <= R[0]) || (x >= X[0]); + else if (X[0] < R[0]) + intersection_outside_XR = (x >= R[0]) || (x <= X[0]); + else if (X[1] > R[1]) + intersection_outside_XR = (y <= R[1]) || (y >= X[1]); + else if (X[1] < R[1]) + intersection_outside_XR = (y >= R[1]) || (y <= X[1]); + else + intersection_outside_XR = true; + } + if (intersection_outside_XR) + return true; + else + { + if (S1[0] > S2[0]) + return (x <= S2[0]) || (x >= S1[0]); + else if (S1[0] < S2[0]) + return (x >= S2[0]) || (x <= S1[0]); + else if (S1[1] > S2[1]) + return (y <= S2[1]) || (y >= S1[1]); + else if (S1[1] < S2[1]) + return (y >= S2[1]) || (y <= S1[1]); + else + return false; + } + } + + /** \brief GreedyProjectionTriangulation is an implementation of a greedy triangulation algorithm for 3D points + * based on local 2D projections. It assumes locally smooth surfaces and relatively smooth transitions between + * areas with different point densities. + * \author Zoltan Csaba Marton + * \ingroup surface + */ + class GreedyProjection : public Filter + { + public: + + enum class GP3Type + { + NONE = -1, // not-defined + FREE = 0, + FRINGE = 1, + BOUNDARY = 2, + COMPLETED = 3 + }; + + /** \brief Empty constructor. */ + GreedyProjection () : + mu_ (0), + search_radius_ (0), // must be set by user + nnn_ (100), + minimum_angle_ (M_PI/18), // 10 degrees + maximum_angle_ (2*M_PI/3), // 120 degrees + eps_angle_(M_PI/4), //45 degrees, + consistent_(false), + consistent_ordering_ (false), + angles_ (), + R_ (), + state_ (), + source_ (), + ffn_ (), + sfn_ (), + part_ (), + fringe_queue_ (), + is_current_free_ (false), + current_index_ (), + prev_is_ffn_ (false), + prev_is_sfn_ (false), + next_is_ffn_ (false), + next_is_sfn_ (false), + changed_1st_fn_ (false), + changed_2nd_fn_ (false), + new2boundary_ (), + already_connected_ (false), + proj_qp_ (), + u_ (), + v_ (), + uvn_ffn_ (), + uvn_sfn_ (), + uvn_next_ffn_ (), + uvn_next_sfn_ (), + tmp_ (), + view_(nullptr) + {}; + + /** \brief Don't consider points for triangulation if their normal deviates more than this value from the query point's normal. + * \param[in] eps_angle maximum surface angle + * \note As normal estimation methods usually give smooth transitions at sharp edges, this ensures correct triangulation + * by avoiding connecting points from one side to points from the other through forcing the use of the edge points. + */ + inline void + setMaximumSurfaceAngle (double eps_angle) { eps_angle_ = eps_angle; } + + /** \brief Get the maximum surface angle. */ + inline double + getMaximumSurfaceAngle () const { return (eps_angle_); } + + /** \brief Set the flag if the input normals are oriented consistently. + * \param[in] consistent set it to true if the normals are + consistently oriented + */ + inline void + setNormalConsistency (bool consistent) { consistent_ = consistent; } + + /** \brief Get the flag for consistently oriented normals. */ + inline bool + getNormalConsistency () const { return (consistent_); } + + /** \brief Set the flag to order the resulting triangle vertices + consistently (positive direction around normal). + * @note Assumes consistently oriented normals (towards the viewpoint) + -- see setNormalConsistency () + * \param[in] consistent_ordering set it to true if triangle + vertices should be ordered consistently + */ + inline void + setConsistentVertexOrdering (bool consistent_ordering) { consistent_ordering_ = consistent_ordering; } + + /** \brief Get the flag signaling consistently ordered triangle vertices. */ + inline bool + getConsistentVertexOrdering () const { return (consistent_ordering_); } + + /** \brief Get the state of each point after reconstruction. + * \note Options are defined as constants: FREE, FRINGE, COMPLETED, BOUNDARY and NONE + */ + /** + inline std::vector + getPointStates () const { return (state_); } + **/ + + /** \brief Get the ID of each point after reconstruction. + * \note parts are numbered from 0, a -1 denotes unconnected points + */ + inline std::vector + getPartIDs () const { return (part_); } + + + /** \brief Get the sfn list. */ + inline std::vector + getSFN () const { return (sfn_); } + + /** \brief Get the ffn list. */ + inline std::vector + getFFN () const { return (ffn_); } + + protected: + /** \brief The nearest neighbor distance multiplier to obtain the + final search radius. + */ + double mu_; + + /** \brief Nearest neighbors search radius for each point and + the maximum edge length. + */ + double search_radius_; + + /** \brief Maximum number of nearest neighbors accepted by searching. */ + int nnn_; + + /** \brief The preferred minimum angle for the triangles. */ + double minimum_angle_; + + /** \brief The maximum angle for the triangles. */ + double maximum_angle_; + + /** \brief Maximum surface angle. */ + double eps_angle_; + + /** \brief Set this to true if the normals of the input are + consistently oriented. + */ + bool consistent_; + + /** \brief Set this to true if the output triangle vertices should + be consistently oriented. + */ + bool consistent_ordering_; + + private: + /** \brief Struct for storing the angles to nearest neighbors **/ + struct nnAngle + { + double angle; + int index; + int nnIndex; + bool visible; + }; + + /** \brief Struct for storing the edges starting from a fringe point **/ + struct doubleEdge + { + doubleEdge () : index (0), first (), second () {} + int index; + Eigen::Vector2f first; + Eigen::Vector2f second; + }; + + struct Triangle + { + PointId a; + PointId b; + PointId c; + }; + + // Variables made global to decrease the number of parameters to helper functions + + /** \brief A list of angles to neighbors **/ + std::vector angles_; + /** \brief Index of the current query point **/ + int R_; + /** \brief List of point states **/ + std::vector state_; + /** \brief List of sources **/ + std::vector source_; + /** \brief List of fringe neighbors in one direction **/ + std::vector ffn_; + /** \brief List of fringe neighbors in other direction **/ + std::vector sfn_; + /** \brief Connected component labels for each point **/ + std::vector part_; + /** \brief Points on the outer edge from which the mesh is grown **/ + std::vector fringe_queue_; + + /** \brief Flag to set if the current point is free **/ + bool is_current_free_; + /** \brief Current point's index **/ + int current_index_; + /** \brief Flag set if the previous point is the first fringe neighbor **/ + bool prev_is_ffn_; + /** \brief Flag to set if the next point is the second fringe neighbor **/ + bool prev_is_sfn_; + /** \brief Flag to set if the next point is the first fringe neighbor **/ + bool next_is_ffn_; + /** \brief Flag to set if the next point is the second fringe neighbor **/ + bool next_is_sfn_; + /** \brief Flag to set if the first fringe neighbor was changed **/ + bool changed_1st_fn_; + /** \brief Flag to set if the second fringe neighbor was changed **/ + bool changed_2nd_fn_; + /** \brief New boundary point **/ + int new2boundary_; + /** \brief Flag to set if the next neighbor was already connected in the previous step. + * To avoid inconsistency it should not be connected again. + */ + bool already_connected_; + + /** \brief Point coordinates projected onto the plane defined by the point normal **/ + Eigen::Vector3d proj_qp_; + /** \brief First coordinate vector of the 2D coordinate frame **/ + Eigen::Vector3d u_; + /** \brief Second coordinate vector of the 2D coordinate frame **/ + Eigen::Vector3d v_; + /** \brief 2D coordinates of the first fringe neighbor **/ + Eigen::Vector2d uvn_ffn_; + /** \brief 2D coordinates of the second fringe neighbor **/ + Eigen::Vector2d uvn_sfn_; + /** \brief 2D coordinates of the 1st fringe neighbor of the next point **/ + Eigen::Vector2d uvn_next_ffn_; + /** \brief 2D coordinates of the 2nd fringe neighbor of the next point **/ + Eigen::Vector2d uvn_next_sfn_; + /** \brief Temporary variable to store 3 coordiantes **/ + Eigen::Vector3d tmp_; + /** \brief Output grid **/ + std::vector grid_; + /** \brief Pointer to current point view. **/ + PointView *view_; + + /** \brief Forms a new triangle by connecting the current neighbor to the query point + * and the previous neighbor + * \param[in] prev_index index of the previous point + * \param[in] next_index index of the next point + * \param[in] next_next_index index of the point after the next one + * \param[in] uvn_current 2D coordinate of the current point + * \param[in] uvn_prev 2D coordinates of the previous point + * \param[in] uvn_next 2D coordinates of the next point + */ + void + connectPoint (const int prev_index, + const int next_index, + const int next_next_index, + const Eigen::Vector2f &uvn_current, + const Eigen::Vector2f &uvn_prev, + const Eigen::Vector2f &uvn_next); + + /** \brief Whenever a query point is part of a boundary loop containing 3 points, that triangle is created + * (called if angle constraints make it possible) + * \param[out] polygons the polygon mesh to be updated + */ + void closeTriangle(); + + /** \brief Get the list of containing triangles for each vertex in a PolygonMesh + * \param[in] polygonMesh the input polygon mesh + */ + /** + std::vector > + getTriangleList (const pcl::PolygonMesh &input); + **/ + + /** \brief Add a new triangle to the current polygon mesh + * \param[in] a index of the first vertex + * \param[in] b index of the second vertex + * \param[in] c index of the third vertex + */ + + /** + inline void + addTriangle (PointId a, PointId b, PointId c) + { + triangle_.vertices.resize (3); + if (consistent_ordering_) + { + const PointInT p = input_->at (indices_->at (a)); + const Eigen::Vector3f pv = p.getVector3fMap (); + if (p.getNormalVector3fMap ().dot ( + (pv - input_->at (indices_->at (b)).getVector3fMap ()).cross ( + pv - input_->at (indices_->at (c)).getVector3fMap ()) ) > 0) + { + triangle_.vertices[0] = a; + triangle_.vertices[1] = b; + triangle_.vertices[2] = c; + } + else + { + triangle_.vertices[0] = a; + triangle_.vertices[1] = c; + triangle_.vertices[2] = b; + } + } + else + { + triangle_.vertices[0] = a; + triangle_.vertices[1] = b; + triangle_.vertices[2] = c; + } + polygons.push_back (triangle_); + } +**/ + + /** \brief Add a new vertex to the advancing edge front and set its source point + * \param[in] v index of the vertex that was connected + * \param[in] s index of the source point + */ + inline void + addFringePoint (int v, int s) + { + source_[v] = s; + part_[v] = part_[s]; + fringe_queue_.push_back(v); + } + + bool stateSet(PointId idx) + { + return state_[idx] != GP3Type::NONE && state_[idx] != GP3Type::FREE; + } + + void addArgs(ProgramArgs& args); + void initialize(); + void filter(PointView& view); + void addTriangle(PointId a, PointId b, PointId c); + Eigen::Vector3d getCoord(PointId id); + Eigen::Vector3d getNormalCoord(PointId id); + }; + +} // namespace pdal + diff --git a/pdal/Dimension.json b/pdal/Dimension.json index ee9c4bd775..a44ab14e25 100644 --- a/pdal/Dimension.json +++ b/pdal/Dimension.json @@ -15,6 +15,21 @@ "description": "Z coordinate" }, { + "name": "NormalX", + "type": "double", + "description": "X component of normal vector" + }, + { + "name": "NormalY", + "type": "double", + "description": "Y component of normal vector" + }, + { + "name": "NormalZ", + "type": "double", + "description": "Z component of normal vector" + }, + { "name": "Intensity", "type": "uint16", "description": "Representation of the pulse return magnitude" From 48761cd7b5d2e8b647738b43ea8c602949558c02 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Wed, 12 Apr 2017 18:01:55 -0500 Subject: [PATCH 04/57] Fix more errors. --- filters/GreedyProjection.cpp | 60 +++++++++++++++++------------------- 1 file changed, 29 insertions(+), 31 deletions(-) diff --git a/filters/GreedyProjection.cpp b/filters/GreedyProjection.cpp index 1e423ebbdb..63b030e57c 100644 --- a/filters/GreedyProjection.cpp +++ b/filters/GreedyProjection.cpp @@ -208,7 +208,7 @@ void GreedyProjection::filter(PointView& view) u_ = nc.cross (v_); // Projecting point onto the surface - const Eigen::Vector3d coord(getCoord(R_)); + Eigen::Vector3d coord(getCoord(R_)); float dist = nc.dot (coord); proj_qp_ = coord - dist * nc; @@ -379,9 +379,9 @@ void GreedyProjection::filter(PointView& view) Eigen::Vector3d coord(getCoord(R_)); // Locating FFN and SFN to adapt distance threshold - double sqr_source_dist = (coord - getCoord[source_[R_]]).squaredNorm (); - double sqr_ffn_dist = (coord - getCoord[ffn_[R_]]).squaredNorm (); - double sqr_sfn_dist = (coord - getCoord[sfn_[R_]]).squaredNorm (); + double sqr_source_dist = (coord - getCoord(source_[R_])).squaredNorm (); + double sqr_ffn_dist = (coord - getCoord(ffn_[R_])).squaredNorm (); + double sqr_sfn_dist = (coord - getCoord(sfn_[R_])).squaredNorm (); double max_sqr_fn_dist = (std::max)(sqr_ffn_dist, sqr_sfn_dist); double sqr_dist_threshold = (std::min)(sqr_max_edge, sqr_mu * sqrDists[1]); if (max_sqr_fn_dist > sqrDists[nnn_-1]) @@ -411,10 +411,10 @@ void GreedyProjection::filter(PointView& view) v_ = nc.unitOrthogonal (); u_ = nc.cross (v_); - const Eigen::Vector3d coord(getCoord(R_)); + const Eigen::Vector3d c(getCoord(R_)); // Projecting point onto the surface - float dist = nc.dot (coord); - proj_qp_ = coord - dist * nc; + float dist = nc.dot (c); + proj_qp_ = c - dist * nc; // Converting coords, calculating angles and saving the projected // near boundary edges @@ -424,7 +424,7 @@ void GreedyProjection::filter(PointView& view) //ABELL - Check this. for (int i = 1; i < nnn_; i++) { - tmp_ = coords_[nnIdx[i]] - proj_qp_; + tmp_ = getCoord(nnIdx[i]) - proj_qp_; uvn_nn[i][0] = tmp_.dot(u_); uvn_nn[i][1] = tmp_.dot(v_); @@ -447,7 +447,7 @@ void GreedyProjection::filter(PointView& view) if ((ffn_[R_] == nnIdx[i]) || (sfn_[R_] == nnIdx[i])) angles_[i].visible = true; bool same_side = true; - const Eigen::Vector3f neighbor_normal = getNormalCoord(nnIdx[i]); + const Eigen::Vector3d neighbor_normal = getNormalCoord(nnIdx[i]); double cosine = nc.dot (neighbor_normal); if (cosine > 1) cosine = 1; if (cosine < -1) cosine = -1; @@ -467,10 +467,10 @@ void GreedyProjection::filter(PointView& view) doubleEdge e; e.index = i; nr_edge++; - tmp_ = getCoord[ffn_[nnIdx[i]]] - proj_qp_; + tmp_ = getCoord(ffn_[nnIdx[i]]) - proj_qp_; e.first[0] = tmp_.dot(u_); e.first[1] = tmp_.dot(v_); - tmp_ = getCoord[sfn_[nnIdx[i]]] - proj_qp_; + tmp_ = getCoord(sfn_[nnIdx[i]]) - proj_qp_; e.second[0] = tmp_.dot(u_); e.second[1] = tmp_.dot(v_); doubleEdges.push_back(e); @@ -563,7 +563,6 @@ void GreedyProjection::filter(PointView& view) } // Sorting angles - auto comp = [](const nnAngle& a1, const nnAngle std::sort(angles_.begin (), angles_.end (), [](const nnAngle& a1, const nnAngle& a2) { @@ -588,13 +587,13 @@ void GreedyProjection::filter(PointView& view) else { if (sqr_max_edge < - (coords_[ffn_[R_]] - coords_[sfn_[R_]]).squaredNorm ()) + (getCoord(ffn_[R_]) - getCoord(sfn_[R_])).squaredNorm ()) { state_[R_] = GP3Type::BOUNDARY; } else { - tmp_ = coords_[source_[R_]] - proj_qp_; + tmp_ = getCoord(source_[R_]) - proj_qp_; uvn_s[0] = tmp_.dot(u_); uvn_s[1] = tmp_.dot(v_); double angleS = atan2(uvn_s[1], uvn_s[0]); @@ -825,7 +824,7 @@ void GreedyProjection::filter(PointView& view) } if ((!gaps[j]) && (sqr_max_edge < - (coords_[angles_[j+1].index] - coords_[angles_[j].index]). + (getCoord(angles_[j+1].index) - getCoord(angles_[j].index)). squaredNorm ())) { gaps[j] = is_boundary = true; @@ -844,8 +843,8 @@ void GreedyProjection::filter(PointView& view) gaps[last_visible] = is_boundary = true; } if ((!gaps[last_visible]) && - (sqr_max_edge < (coords_[angles_[0].index] - - coords_[angles_[last_visible].index]).squaredNorm ())) + (sqr_max_edge < (getCoord(angles_[0].index) - + getCoord(angles_[last_visible].index)).squaredNorm ())) { gaps[last_visible] = is_boundary = true; } @@ -863,8 +862,8 @@ void GreedyProjection::filter(PointView& view) gaps[j] = is_boundary = true; } if ((!gaps[j]) && - (sqr_max_edge < (coords_[angles_[j+1].index] - - coords_[angles_[j].index]).squaredNorm ())) + (sqr_max_edge < (getCoord(angles_[j+1].index) - + getCoord(angles_[j].index)).squaredNorm ())) { gaps[j] = is_boundary = true; } @@ -887,8 +886,8 @@ void GreedyProjection::filter(PointView& view) gaps[j] = is_boundary = true; } if ((!gaps[j]) && - (sqr_max_edge < (coords_[angles_[j+1].index] - - coords_[angles_[j].index]).squaredNorm ())) + (sqr_max_edge < (getCoord(angles_[j+1].index) - + getCoord(angles_[j].index)).squaredNorm ())) { gaps[j] = is_boundary = true; } @@ -962,8 +961,8 @@ void GreedyProjection::filter(PointView& view) else { int erased_idx = static_cast (to_erase.size ()) -1; - for (auto prev_it = it-1; - (erased_idx != -1) && (it != angleIdx.begin()); it--) + auto prev_it = it - 1; + for (; (erased_idx != -1) && (it != angleIdx.begin()); it--) if (*it == to_erase[erased_idx]) erased_idx--; else @@ -1175,13 +1174,13 @@ void GreedyProjection::filter(PointView& view) } } log()->get(LogLevel::Debug) << "Number of triangles: " << grid_.size() << - ".\n"); + ".\n"; log()->get(LogLevel::Debug) << "Number of unconnected parts: " << nr_parts << - ".\n"); + ".\n"; if (increase_nnn4fn > 0) log()->get(LogLevel::Warning) << "Number of neighborhood size " "increase requests for fringe neighbors: " << increase_nnn4fn << - ".\n"); + ".\n"; if (increase_nnn4s > 0) log()->get(LogLevel::Warning) << "Number of neighborhood size " "increase requests for source: " << increase_nnn4s << ".\n"; @@ -1196,7 +1195,6 @@ void GreedyProjection::filter(PointView& view) log()->get(LogLevel::Debug) << "Number of processed points: " << fringe_queue_.size() << " / " << view.size() << "!\n"; view_ = nullptr; - return (true); } GreedyProjection::closeTriangle () @@ -1590,10 +1588,10 @@ void GreedyProjection::connectPoint ( } else { - tmp_ = coords_[ffn_[next_index]] - proj_qp_; + tmp_ = getCoord(ffn_[next_index]) - proj_qp_; uvn_next_ffn_[0] = tmp_.dot(u_); uvn_next_ffn_[1] = tmp_.dot(v_); - tmp_ = coords_[sfn_[next_index]] - proj_qp_; + tmp_ = getCoord(sfn_[next_index]) - proj_qp_; uvn_next_sfn_[0] = tmp_.dot(u_); uvn_next_sfn_[1] = tmp_.dot(v_); @@ -1740,10 +1738,10 @@ void GreedyProjection::connectPoint ( } else { - tmp_ = coords_[ffn_[next_index]] - proj_qp_; + tmp_ = getCoord(ffn_[next_index]) - proj_qp_; uvn_next_ffn_[0] = tmp_.dot(u_); uvn_next_ffn_[1] = tmp_.dot(v_); - tmp_ = coords_[sfn_[next_index]] - proj_qp_; + tmp_ = getCoord(sfn_[next_index]) - proj_qp_; uvn_next_sfn_[0] = tmp_.dot(u_); uvn_next_sfn_[1] = tmp_.dot(v_); From 727f94f3f7e00dbd9f053c9bca4ee1be6d2faca1 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Thu, 13 Apr 2017 12:39:32 -0500 Subject: [PATCH 05/57] It compiles. --- filters/GreedyProjection.cpp | 149 +++++++++++++++++------------------ filters/GreedyProjection.hpp | 64 ++++++--------- 2 files changed, 97 insertions(+), 116 deletions(-) diff --git a/filters/GreedyProjection.cpp b/filters/GreedyProjection.cpp index 63b030e57c..b56ca9a540 100644 --- a/filters/GreedyProjection.cpp +++ b/filters/GreedyProjection.cpp @@ -118,10 +118,10 @@ void GreedyProjection::filter(PointView& view) int part_index = 0; // 2D coordinates of points - const Eigen::Vector2f uvn_nn_qp_zero = Eigen::Vector2f::Zero(); - Eigen::Vector2f uvn_current; - Eigen::Vector2f uvn_prev; - Eigen::Vector2f uvn_next; + const Eigen::Vector2d uvn_nn_qp_zero = Eigen::Vector2d::Zero(); + Eigen::Vector2d uvn_current; + Eigen::Vector2d uvn_prev; + Eigen::Vector2d uvn_next; // initializing fields already_connected_ = false; // see declaration for comments :P @@ -132,11 +132,11 @@ void GreedyProjection::filter(PointView& view) source_.clear (); ffn_.clear (); sfn_.clear (); - part_.resize(view.size (), -1); // indices of point's part + part_.resize(view.size ()); // indices of point's part state_.resize(view.size (), GP3Type::FREE); - source_.resize(view.size (), -1); - ffn_.resize(view.size (), -1); - sfn_.resize(view.size (), -1); + source_.resize(view.size ()); + ffn_.resize(view.size()); + sfn_.resize(view.size()); fringe_queue_.clear (); int fqIdx = 0; // current fringe's index in the queue to be processed @@ -172,8 +172,8 @@ void GreedyProjection::filter(PointView& view) int is_free=0, nr_parts=0, increase_nnn4fn=0, increase_nnn4s=0, increase_dist=0, nr_touched = 0; bool is_fringe; angles_.resize(nnn_); - std::vector > uvn_nn (nnn_); - Eigen::Vector2f uvn_s; + std::vector > uvn_nn (nnn_); + Eigen::Vector2d uvn_s; // iterating through fringe points and finishing them until everything is done while (!done) @@ -214,7 +214,7 @@ void GreedyProjection::filter(PointView& view) // Converting coords, calculating angles and saving the // projected near boundary edges - int nr_edge = 0; + size_t nr_edge = 0; std::vector doubleEdges; // nearest neighbor with index 0 is the query point R_ itself //ABELL - Need to check that the above is true. @@ -229,13 +229,10 @@ void GreedyProjection::filter(PointView& view) angles_[i].angle = atan2(uvn_nn[i][1], uvn_nn[i][0]); // initializing angle descriptors angles_[i].index = nnIdx[i]; - if ( - (state_[nnIdx[i]] == GP3Type::COMPLETED) || + if ((state_[nnIdx[i]] == GP3Type::COMPLETED) || (state_[nnIdx[i]] == GP3Type::BOUNDARY) || (state_[nnIdx[i]] == GP3Type::NONE) || - (nnIdx[i] == -1) /// NOTE: discarding NaN points and those that are not in indices_ - || (sqrDists[i] > sqr_dist_threshold) - ) + (sqrDists[i] > sqr_dist_threshold)) angles_[i].visible = false; else angles_[i].visible = true; @@ -265,16 +262,16 @@ void GreedyProjection::filter(PointView& view) (sfn_[R_] != nnIdx[i])) { bool visibility = true; - for (int j = 0; j < nr_edge; j++) + for (size_t j = 0; j < nr_edge; j++) { if (ffn_[nnIdx[doubleEdges[j].index]] != nnIdx[i]) visibility = isVisible(uvn_nn[i], uvn_nn[doubleEdges[j].index], - doubleEdges[j].first, Eigen::Vector2f::Zero()); + doubleEdges[j].first, Eigen::Vector2d::Zero()); if (!visibility) break; if (sfn_[nnIdx[doubleEdges[j].index]] != nnIdx[i]) visibility = isVisible(uvn_nn[i], uvn_nn[doubleEdges[j].index], - doubleEdges[j].second, Eigen::Vector2f::Zero()); + doubleEdges[j].second, Eigen::Vector2d::Zero()); if (!visibility == false) break; } @@ -286,10 +283,10 @@ void GreedyProjection::filter(PointView& view) PointId left = 1; do { - while ((left < nnn_) && + while ((left < (PointId)nnn_) && ((!angles_[left].visible) || stateSet(nnIdx[left]))) left++; - if (left >= nnn_) + if (left >= (PointId)nnn_) break; else { @@ -418,7 +415,7 @@ void GreedyProjection::filter(PointView& view) // Converting coords, calculating angles and saving the projected // near boundary edges - int nr_edge = 0; + size_t nr_edge = 0; std::vector doubleEdges; // nearest neighbor with index 0 is the query point R_ itself //ABELL - Check this. @@ -434,13 +431,10 @@ void GreedyProjection::filter(PointView& view) // initializing angle descriptors angles_[i].index = nnIdx[i]; angles_[i].nnIndex = i; - if ( - (state_[nnIdx[i]] == GP3Type::COMPLETED) || + if ((state_[nnIdx[i]] == GP3Type::COMPLETED) || (state_[nnIdx[i]] == GP3Type::BOUNDARY) || (state_[nnIdx[i]] == GP3Type::NONE) || - (nnIdx[i] == -1) /// NOTE: discarding NaN points and those that are not in indices_ - || (sqrDists[i] > sqr_dist_threshold) - ) + (sqrDists[i] > sqr_dist_threshold)) angles_[i].visible = false; else angles_[i].visible = true; @@ -540,21 +534,24 @@ void GreedyProjection::filter(PointView& view) (ffn_[R_] != nnIdx[i]) && (sfn_[R_] != nnIdx[i])) { bool visibility = true; - for (int j = 0; j < nr_edge; j++) + for (size_t j = 0; j < nr_edge; j++) { - if (doubleEdges[j].index != i) + //ABELL - This seems weird. i is just a count of the nearest + // neighbors, not a point index. Are some indicies the index + // of the nearest neighbors? If so, confusing. + if (doubleEdges[j].index != (PointId)i) { - int f = ffn_[nnIdx[doubleEdges[j].index]]; + PointId f = ffn_[nnIdx[doubleEdges[j].index]]; if ((f != nnIdx[i]) && (f != R_)) visibility = isVisible(uvn_nn[i], uvn_nn[doubleEdges[j].index], - doubleEdges[j].first, Eigen::Vector2f::Zero()); + doubleEdges[j].first, Eigen::Vector2d::Zero()); if (visibility == false) break; - int s = sfn_[nnIdx[doubleEdges[j].index]]; + PointId s = sfn_[nnIdx[doubleEdges[j].index]]; if ((s != nnIdx[i]) && (s != R_)) visibility = isVisible(uvn_nn[i], uvn_nn[doubleEdges[j].index], - doubleEdges[j].second, Eigen::Vector2f::Zero()); + doubleEdges[j].second, Eigen::Vector2d::Zero()); if (visibility == false) break; } @@ -923,9 +920,9 @@ void GreedyProjection::filter(PointView& view) double angle_so_far = 0, angle_would_be; double max_combined_angle = (std::min)(maximum_angle_, M_PI-2*minimum_angle_); - Eigen::Vector2f X; - Eigen::Vector2f S1; - Eigen::Vector2f S2; + Eigen::Vector2d X; + Eigen::Vector2d S1; + Eigen::Vector2d S2; std::vector to_erase; for (auto it = angleIdx.begin()+1; it != angleIdx.end()-1; it++) { @@ -1197,7 +1194,7 @@ void GreedyProjection::filter(PointView& view) view_ = nullptr; } -GreedyProjection::closeTriangle () +void GreedyProjection::closeTriangle () { state_[R_] = GP3Type::COMPLETED; addTriangle (angles_[0].index, angles_[1].index, R_); @@ -1229,7 +1226,7 @@ GreedyProjection::closeTriangle () } void GreedyProjection::connectPoint ( - const int prev_index, const int next_index, const int next_next_index, + PointId prev_index, PointId next_index, PointId next_next_index, const Eigen::Vector2d &uvn_current, const Eigen::Vector2d &uvn_prev, const Eigen::Vector2d &uvn_next) @@ -1335,14 +1332,14 @@ void GreedyProjection::connectPoint ( if (prev_ffn && next_sfn && prev_sfn && next_ffn) { /* should be never the case */ - double prev2f = (getCoord[ffn_[current_index_]] - - getCoord[prev_index]).squaredNorm (); - double next2s = (getCoord[sfn_[current_index_]] - - getCoord[next_index]).squaredNorm (); - double prev2s = (getCoord[sfn_[current_index_]] - - getCoord[prev_index]).squaredNorm (); - double next2f = (getCoord[ffn_[current_index_]] - - getCoord[next_index]).squaredNorm (); + double prev2f = (getCoord(ffn_[current_index_]) - + getCoord(prev_index)).squaredNorm (); + double next2s = (getCoord(sfn_[current_index_]) - + getCoord(next_index)).squaredNorm (); + double prev2s = (getCoord(sfn_[current_index_]) - + getCoord(prev_index)).squaredNorm (); + double next2f = (getCoord(ffn_[current_index_]) - + getCoord(next_index)).squaredNorm (); if (prev2f < prev2s) { if (prev2f < next2f) @@ -1381,10 +1378,10 @@ void GreedyProjection::connectPoint ( else if (prev_ffn && next_sfn) { /* a clear case */ - double prev2f = (getCoord[ffn_[current_index_]] - - getCoord[prev_index]).squaredNorm (); - double next2s = (getCoord[sfn_[current_index_]] - - getCoord[next_index]).squaredNorm (); + double prev2f = (getCoord(ffn_[current_index_]) - + getCoord(prev_index)).squaredNorm(); + double next2s = (getCoord(sfn_[current_index_]) - + getCoord(next_index)).squaredNorm (); if (prev2f < next2s) min_dist = 0; else @@ -1393,10 +1390,10 @@ void GreedyProjection::connectPoint ( else if (prev_sfn && next_ffn) { /* a clear case */ - double prev2s = (getCoord[sfn_[current_index_]] - - getCoord[prev_index]).squaredNorm (); - double next2f = (getCoord[ffn_[current_index_]] - - getCoord[next_index]).squaredNorm (); + double prev2s = (getCoord(sfn_[current_index_]) - + getCoord(prev_index)).squaredNorm (); + double next2f = (getCoord(ffn_[current_index_]) - + getCoord(next_index)).squaredNorm (); if (prev2s < next2f) min_dist = 1; else @@ -1414,12 +1411,12 @@ void GreedyProjection::connectPoint ( /* messed up cases */ else if (prev_ffn) { - double prev2f = (getCoord[ffn_[current_index_]] - - getCoord[prev_index]).squaredNorm (); + double prev2f = (getCoord(ffn_[current_index_]) - + getCoord(prev_index)).squaredNorm (); if (prev_sfn) { - double prev2s = (getCoord[sfn_[current_index_]] - - getCoord[prev_index]).squaredNorm (); + double prev2s = (getCoord(sfn_[current_index_]) - + getCoord(prev_index)).squaredNorm (); if (prev2s < prev2f) min_dist = 1; else @@ -1427,8 +1424,8 @@ void GreedyProjection::connectPoint ( } else if (next_ffn) { - double next2f = (getCoord[ffn_[current_index_]] - - getCoord[next_index]).squaredNorm (); + double next2f = (getCoord(ffn_[current_index_]) - + getCoord(next_index)).squaredNorm (); if (next2f < prev2f) min_dist = 2; else @@ -1437,12 +1434,12 @@ void GreedyProjection::connectPoint ( } else if (next_sfn) { - double next2s = (getCoord[sfn_[current_index_]] - - getCoord[next_index]).squaredNorm (); + double next2s = (getCoord(sfn_[current_index_]) - + getCoord(next_index)).squaredNorm (); if (prev_sfn) { - double prev2s = (getCoord[sfn_[current_index_]] - - getCoord[prev_index]).squaredNorm (); + double prev2s = (getCoord(sfn_[current_index_]) - + getCoord(prev_index)).squaredNorm (); if (prev2s < next2s) min_dist = 1; else @@ -1450,8 +1447,8 @@ void GreedyProjection::connectPoint ( } else if (next_ffn) { - double next2f = (getCoord[ffn_[current_index_]] - - getCoord[next_index]).squaredNorm (); + double next2f = (getCoord(ffn_[current_index_]) - + getCoord(next_index)).squaredNorm (); if (next2f < next2s) min_dist = 2; else @@ -1549,7 +1546,7 @@ void GreedyProjection::connectPoint ( case 2://next2f: { addTriangle (current_index_, ffn_[current_index_], next_index); - int neighbor_update = next_index; + PointId neighbor_update = next_index; /* updating next_index */ if (!stateSet(next_index)) @@ -1605,10 +1602,10 @@ void GreedyProjection::connectPoint ( int connect2ffn = -1; if (ffn_next_ffn && sfn_next_ffn) { - double fn2f = (getCoord[ffn_[current_index_]] - - getCoord[ffn_[next_index]]).squaredNorm (); - double sn2f = (getCoord[ffn_[current_index_]] - - getCoord[sfn_[next_index]]).squaredNorm (); + double fn2f = (getCoord(ffn_[current_index_]) - + getCoord(ffn_[next_index])).squaredNorm (); + double sn2f = (getCoord(ffn_[current_index_]) - + getCoord(sfn_[next_index])).squaredNorm (); if (fn2f < sn2f) connect2ffn = 0; else @@ -1699,7 +1696,7 @@ void GreedyProjection::connectPoint ( case 3://next2s: { addTriangle (current_index_, sfn_[current_index_], next_index); - int neighbor_update = next_index; + PointId neighbor_update = next_index; /* updating next_index */ if (!stateSet(next_index)) @@ -1755,10 +1752,10 @@ void GreedyProjection::connectPoint ( int connect2sfn = -1; if (ffn_next_sfn && sfn_next_sfn) { - double fn2s = (getCoord[sfn_[current_index_]] - - getCoord[ffn_[next_index]]).squaredNorm (); - double sn2s = (getCoord[sfn_[current_index_]] - - getCoord[sfn_[next_index]]).squaredNorm (); + double fn2s = (getCoord(sfn_[current_index_]) - + getCoord(ffn_[next_index])).squaredNorm (); + double sn2s = (getCoord(sfn_[current_index_]) - + getCoord(sfn_[next_index])).squaredNorm (); if (fn2s < sn2s) connect2sfn = 0; else diff --git a/filters/GreedyProjection.hpp b/filters/GreedyProjection.hpp index 79dc8a2fad..66e1e59f4f 100644 --- a/filters/GreedyProjection.hpp +++ b/filters/GreedyProjection.hpp @@ -57,7 +57,7 @@ namespace pdal inline bool isVisible (const Eigen::Vector2d &X, const Eigen::Vector2d &S1, const Eigen::Vector2d &S2, - const Eigen::Vector2d &R = Eigen::Vector2f::Zero ()) + const Eigen::Vector2d &R = Eigen::Vector2d::Zero ()) { double a0 = S1[1] - S2[1]; double b0 = S2[0] - S1[0]; @@ -65,7 +65,7 @@ namespace pdal double a1 = -X[1]; double b1 = X[0]; double c1 = 0; - if (R != Eigen::Vector2f::Zero()) + if (R != Eigen::Vector2d::Zero()) { a1 += R[1]; b1 -= R[0]; @@ -76,7 +76,7 @@ namespace pdal double y = (a1*c0 - a0*c1) / div; bool intersection_outside_XR; - if (R == Eigen::Vector2f::Zero()) + if (R == Eigen::Vector2d::Zero()) { if (X[0] > 0) intersection_outside_XR = (x <= 0) || (x >= X[0]); @@ -217,25 +217,6 @@ namespace pdal /** \brief Get the state of each point after reconstruction. * \note Options are defined as constants: FREE, FRINGE, COMPLETED, BOUNDARY and NONE */ - /** - inline std::vector - getPointStates () const { return (state_); } - **/ - - /** \brief Get the ID of each point after reconstruction. - * \note parts are numbered from 0, a -1 denotes unconnected points - */ - inline std::vector - getPartIDs () const { return (part_); } - - - /** \brief Get the sfn list. */ - inline std::vector - getSFN () const { return (sfn_); } - - /** \brief Get the ffn list. */ - inline std::vector - getFFN () const { return (ffn_); } protected: /** \brief The nearest neighbor distance multiplier to obtain the @@ -275,7 +256,7 @@ namespace pdal struct nnAngle { double angle; - int index; + PointId index; int nnIndex; bool visible; }; @@ -284,13 +265,16 @@ namespace pdal struct doubleEdge { doubleEdge () : index (0), first (), second () {} - int index; - Eigen::Vector2f first; - Eigen::Vector2f second; + PointId index; + Eigen::Vector2d first; + Eigen::Vector2d second; }; struct Triangle { + Triangle(PointId ta, PointId tb, PointId tc) : a(ta), b(tb), c(tc) + {} + PointId a; PointId b; PointId c; @@ -301,24 +285,24 @@ namespace pdal /** \brief A list of angles to neighbors **/ std::vector angles_; /** \brief Index of the current query point **/ - int R_; + PointId R_; /** \brief List of point states **/ std::vector state_; /** \brief List of sources **/ - std::vector source_; + std::vector source_; /** \brief List of fringe neighbors in one direction **/ - std::vector ffn_; + std::vector ffn_; /** \brief List of fringe neighbors in other direction **/ - std::vector sfn_; + std::vector sfn_; /** \brief Connected component labels for each point **/ - std::vector part_; + std::vector part_; /** \brief Points on the outer edge from which the mesh is grown **/ - std::vector fringe_queue_; + std::vector fringe_queue_; /** \brief Flag to set if the current point is free **/ bool is_current_free_; /** \brief Current point's index **/ - int current_index_; + PointId current_index_; /** \brief Flag set if the previous point is the first fringe neighbor **/ bool prev_is_ffn_; /** \brief Flag to set if the next point is the second fringe neighbor **/ @@ -332,7 +316,7 @@ namespace pdal /** \brief Flag to set if the second fringe neighbor was changed **/ bool changed_2nd_fn_; /** \brief New boundary point **/ - int new2boundary_; + PointId new2boundary_; /** \brief Flag to set if the next neighbor was already connected in the previous step. * To avoid inconsistency it should not be connected again. */ @@ -369,12 +353,12 @@ namespace pdal * \param[in] uvn_next 2D coordinates of the next point */ void - connectPoint (const int prev_index, - const int next_index, - const int next_next_index, - const Eigen::Vector2f &uvn_current, - const Eigen::Vector2f &uvn_prev, - const Eigen::Vector2f &uvn_next); + connectPoint (PointId prev_index, + PointId next_index, + PointId next_next_index, + const Eigen::Vector2d &uvn_current, + const Eigen::Vector2d &uvn_prev, + const Eigen::Vector2d &uvn_next); /** \brief Whenever a query point is part of a boundary loop containing 3 points, that triangle is created * (called if angle constraints make it possible) From 2b1819473fb5ac433284c2d89b6090a0f8041cbf Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Thu, 13 Apr 2017 14:07:29 -0500 Subject: [PATCH 06/57] Add bounds option to writers.gdal doc. --- doc/stages/writers.gdal.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/doc/stages/writers.gdal.rst b/doc/stages/writers.gdal.rst index 0ff2baa4b3..18efa48174 100644 --- a/doc/stages/writers.gdal.rst +++ b/doc/stages/writers.gdal.rst @@ -135,3 +135,12 @@ window_size dimension A dimension name to use for the interpolation. [Default: ``Z``] + +.. _bounds: + +bounds + The bounds of the data to be written. Points not in bounds are discarded. + The format is ([minx, maxx],[miny,maxy]). + +.. note:: + The bounds_ option is required when a pipeline is run in streaming mode. From 57d0c3da0e1d2d1fec069067124ea5995080d77d Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Fri, 14 Apr 2017 07:56:50 -0500 Subject: [PATCH 07/57] Remove code moved to PointView.cpp. --- pdal/PointView.hpp | 27 +-------------------------- 1 file changed, 1 insertion(+), 26 deletions(-) diff --git a/pdal/PointView.hpp b/pdal/PointView.hpp index 3031b0b014..75687c7696 100644 --- a/pdal/PointView.hpp +++ b/pdal/PointView.hpp @@ -95,6 +95,7 @@ class PDAL_DLL PointView : public PointContainer { // We use size() instead of the index end because temp points // might have been placed at the end of the buffer. + // We're essentially ditching temp points. auto thisEnd = m_index.begin() + size(); auto bufEnd = buf.m_index.begin() + buf.size(); m_index.insert(thisEnd, buf.m_index.begin(), bufEnd); @@ -528,32 +529,6 @@ void PointView::setField(Dimension::Id dim, PointId idx, T val) } } -/** -void PointView::setFieldInternal(Dimension::Id dim, PointId idx, - const void *value) -{ - PointId rawId = 0; - if (idx == size()) - { - rawId = m_pointTable.addPoint(); - m_index.push_back(rawId); - m_size++; - assert(m_temps.empty()); - } - else if (idx > size()) - { - std::cerr << "Point index must increment.\n"; - //error - throw? - return; - } - else - { - rawId = m_index[idx]; - } - m_pointTable.setFieldInternal(dim, rawId, value); -} -**/ - inline void PointView::appendPoint(const PointView& buffer, PointId id) { // Invalid 'id' is a programmer error. From 85c05712b8b07dc90b7090ac1374bb747720a6b0 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Fri, 14 Apr 2017 10:52:12 -0500 Subject: [PATCH 08/57] Allow longname synonyms in ProgramArgs. Close #1560 --- pdal/util/ProgramArgs.hpp | 19 +++++++++++++++++++ plugins/sqlite/io/SQLiteWriter.cpp | 1 + test/unit/ProgramArgsTest.cpp | 17 +++++++++++++++++ 3 files changed, 37 insertions(+) diff --git a/pdal/util/ProgramArgs.hpp b/pdal/util/ProgramArgs.hpp index cf4dc49c6b..19542b0751 100644 --- a/pdal/util/ProgramArgs.hpp +++ b/pdal/util/ProgramArgs.hpp @@ -1159,6 +1159,24 @@ class ProgramArgs } } + /** + Add a synonym for an argument. + + \param name Longname of existing arugment. + \param synonym Synonym for argument. + */ + void addSynonym(const std::string& name, const std::string& synonym) + { + Arg *arg = findLongArg(name); + if (!arg) + throw arg_error("Can't set synonym for argument '" + name + "'. " + "Argument doesn't exist."); + if (synonym.empty()) + throw arg_error("Invalid (empty) synonym for argument '" + + name + "'."); + addLongArg(synonym, arg); + } + /** Reset the state of all arguments and bound variables as if no parsing had occurred. @@ -1312,6 +1330,7 @@ class ProgramArgs } out << "]"; } + private: /* Split an argument name into longname and shortname. diff --git a/plugins/sqlite/io/SQLiteWriter.cpp b/plugins/sqlite/io/SQLiteWriter.cpp index 46fb587b6f..4c005dc0af 100644 --- a/plugins/sqlite/io/SQLiteWriter.cpp +++ b/plugins/sqlite/io/SQLiteWriter.cpp @@ -75,6 +75,7 @@ void SQLiteWriter::addArgs(ProgramArgs& args) m_cloud_table).setPositional(); args.add("connection", "SQL connection string", m_connection).setPositional(); + args.addSynonym("connection", "filename"); args.add("cloud_column_name", "Cloud column name", m_cloud_column, "id"); args.add("module", "Module name", m_modulename); args.add("srid", "SRID", m_srid, 4326U); diff --git a/test/unit/ProgramArgsTest.cpp b/test/unit/ProgramArgsTest.cpp index 6ca8fbbd47..4800db66ec 100644 --- a/test/unit/ProgramArgsTest.cpp +++ b/test/unit/ProgramArgsTest.cpp @@ -228,6 +228,23 @@ TEST(ProgramArgsTest, t4) HANDLE_EXCEPTION(args.add("", "Foo description", m_foo, "foo")); } +TEST(ProgramArgsTest, synonym) +{ + ProgramArgs args; + + std::string m_foo; + + args.add("foo,f", "Foo description", m_foo, "foo"); + args.addSynonym("foo", "bar"); + StringList s = toStringList("--bar"); + EXPECT_THROW(args.parse(s), arg_error); + + s = toStringList("--bar=TestFoo"); + args.reset(); + args.parse(s); + EXPECT_EQ(m_foo, "TestFoo"); +} + TEST(ProgramArgsTest, positional) { ProgramArgs args; From 22bae5ea114b4354314add4643fa656e767e9448 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Fri, 14 Apr 2017 16:17:52 -0500 Subject: [PATCH 09/57] Make sure cell in which point itself lies is in bounds. Close #1563 --- io/GDALGrid.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/io/GDALGrid.cpp b/io/GDALGrid.cpp index af38bcd89e..4511c02d3b 100644 --- a/io/GDALGrid.cpp +++ b/io/GDALGrid.cpp @@ -288,7 +288,9 @@ void GDALGrid::addPoint(double x, double y, double z) // This is a questionable case. If a point is in a cell, shouldn't // it just be counted? double d = distance(iOrigin, jOrigin, x, y); - if (d < m_radius) + if (d < m_radius && + iOrigin >= 0 && jOrigin >= 0 && + iOrigin < (int)m_width && jOrigin <= (int)m_height) update(iOrigin, jOrigin, z, d); } From 8f9c895b47ae8525ddaa779f4054510b8b171e5b Mon Sep 17 00:00:00 2001 From: Bradley J Chambers Date: Mon, 17 Apr 2017 09:57:15 -0400 Subject: [PATCH 10/57] Resolve some dangling references to deprecated extract and classify options extract and classify have been used in several filters, namely outlier and PMF. With the release of PDAL v1.5, we removed these options, opting to always classify points, and to force users to ignore or extract them as they see fit in their pipelines. The outlier filter's extract option can be replicated by adding a range filter stage as shown following the outlier filter: { "type":"filters.range", "limits":"Classification![7:7]" } The PMF filter's extract option can be replicated by adding the following range filter: { "type":"filters.range", "limits":"Classification[2:2]" } The intention here was to remove any ambiguity as to how individual filters implement "extract", and to instead be explicit in the pipeline specification as to how points are extracted. --- doc/apps/translate.rst | 12 ++--- doc/pipeline.rst | 14 +++--- doc/stages/filters.outlier.rst | 29 +++++------ doc/stages/filters.pmf.rst | 50 ++++++++++++------- .../exercises/analysis/denoising/denoise.json | 3 +- .../analysis/denoising/denoising.rst | 21 +++++--- .../ground/ground-run-ground-only.txt | 2 +- .../exercises/analysis/ground/ground.rst | 17 ++++--- .../ground/translate-run-ground-only.txt | 5 +- 9 files changed, 90 insertions(+), 63 deletions(-) diff --git a/doc/apps/translate.rst b/doc/apps/translate.rst index 281d83cf00..7b32bb29a4 100644 --- a/doc/apps/translate.rst +++ b/doc/apps/translate.rst @@ -80,15 +80,16 @@ Example 2: -------------------------------------------------------------------------------- Given these tools, we can now construct a custom pipeline on-the-fly. The -example below uses a simple LAS reader and writer, but stages a PCL-based -voxel grid filter, followed by the PMF filter. We can even set +example below uses a simple LAS reader and writer, but stages a PCL-based voxel +grid filter, followed by the PMF filter and a range filter. We can even set stage-specific parameters as shown. :: - $ pdal translate input.las output.las pclblock pmf \ - --filters.pclblock.json="{\"pipeline\":{\"filters\":[{\"name\":\"VoxelGrid\"}]}}" \ - --filters.pmf.approximate=true --filters.pmf.extract=true + $ pdal translate input.las output.las pclblock pmf range \ + --filters.pclblock.methods="[{\"name\":\"VoxelGrid\"}]" \ + --filters.pmf.approximate=true \ + --filters.range.limits="Classification[2:2]" Example 3: -------------------------------------------------------------------------------- @@ -113,4 +114,3 @@ reference system and writes the result to the file "output.las". $ pdal translate input.las output.las -f filters.reprojection \ --filters.reprojection.out_srs="EPSG:4326" - diff --git a/doc/pipeline.rst b/doc/pipeline.rst index 7a85d99899..97fcedbf06 100644 --- a/doc/pipeline.rst +++ b/doc/pipeline.rst @@ -252,8 +252,9 @@ DTM A common task is to create a digital terrain model (DTM) from the input point cloud. This pipeline infers the reader type, applies an approximate ground -segmentation filter using :ref:`filters.pmf`, and then creates the DTM using -the :ref:`writers.gdal` with only the ground returns. +segmentation filter using :ref:`filters.pmf`, filters out all points but the +ground returns (classification value of 2) using the :ref:`filters.range`, and +then creates the DTM using the :ref:`writers.gdal`. .. code-block:: json @@ -267,9 +268,11 @@ the :ref:`writers.gdal` with only the ground returns. "slope":1.0, "max_distance":2.5, "initial_distance":0.15, - "cell_size":1.0, - "extract":true, - "classify":false + "cell_size":1.0 + }, + { + "type":"range", + "limits":"Classification[2:2]" }, { "type":"writers.gdal", @@ -460,4 +463,3 @@ PDAL. Readers follow the pattern of :ref:`readers.las` or Issuing the command ``pdal info --options`` will list all available stages and their options. See :ref:`info_command` for more. - diff --git a/doc/stages/filters.outlier.rst b/doc/stages/filters.outlier.rst index 3367ad016e..932dc17d65 100644 --- a/doc/stages/filters.outlier.rst +++ b/doc/stages/filters.outlier.rst @@ -7,6 +7,18 @@ filters.outlier The Outlier filter provides two outlier filtering methods: radius and statistical. These two approaches are discussed in further detail below. +It is worth noting that both filtering methods simply apply a classification +value of 7 to the noise points (per the LAS specification). To remove the noise +points altogether, users can add a :ref:`range filter` to their +pipeline, downstream from the outlier filter. + +.. code-block:: json + + { + "type":"filters.range", + "limits":"Classification![7:7]" + } + Statistical Method ------------------------------------------------------------------------------- @@ -39,10 +51,6 @@ We now interate over the pre-computed mean distances :math:`\mu_i` and compare t \text{false,} \phantom{true,} \text{otherwise} \\ \end{cases} -The ``classify`` and ``extract`` options are used to control whether outlier -points are labeled as noise, or removed from the output ``PointView`` -completely. - .. figure:: filters.statisticaloutlier.img1.png :scale: 70 % :alt: Points before outlier removal @@ -96,10 +104,6 @@ of neighbors specified by ``min_k``, it is marked as an outlier. \text{false,} \phantom{true,} \text{otherwise} \\ \end{cases} -The ``classify`` and ``extract`` options are used to control whether outlier -points are labeled as noise, or removed from the output ``PointView`` -completely. - Example ............................................................................... @@ -124,6 +128,9 @@ four neighbors within a radius of 1.0. Options ------------------------------------------------------------------------------- +class + The classification value to apply to outliers. [Default: **7**] + method The outlier removal method. [Default: **statistical**] @@ -138,9 +145,3 @@ mean_k multiplier Standard deviation threshold (statistical method only). [Default: **2.0**] - -classify - Apply classification value of 18 (LAS high noise)? [Default: **true**] - -extract - Extract inlier returns only? [Default: **false**] diff --git a/doc/stages/filters.pmf.rst b/doc/stages/filters.pmf.rst index ab93fe883f..6046843303 100644 --- a/doc/stages/filters.pmf.rst +++ b/doc/stages/filters.pmf.rst @@ -37,25 +37,36 @@ Notes ``initial_distance`` large enough to not exclude these points from the ground. * For a given iteration, the height threshold is determined by multiplying - ``slope`` by ``cell_size`` by the difference in window size between the current - and last iteration, plus the ``initial_distance``. This height threshold is - constant across all cells and is maxed out at the ``max_distance`` value. If - the difference in elevation between a point and its “opened†value (from the - morphological operator) exceeds the height threshold, it is treated as - non-ground. So, bigger slope leads to bigger height thresholds, and these - grow with each iteration (not to exceed the max). With flat terrain, - keep this low, the thresholds are small, and stuff is more aggressively - dumped into non-ground class. In rugged terrain, open things up + ``slope`` by ``cell_size`` by the difference in window size between the + current and last iteration, plus the ``initial_distance``. This height + threshold is constant across all cells and is maxed out at the + ``max_distance`` value. If the difference in elevation between a point and its + “opened†value (from the morphological operator) exceeds the height threshold, + it is treated as non-ground. So, bigger slope leads to bigger height + thresholds, and these grow with each iteration (not to exceed the max). With + flat terrain, keep this low, the thresholds are small, and stuff is more + aggressively dumped into non-ground class. In rugged terrain, open things up a little, but then you can start missing buildings, veg, etc. * Very large ``max_window_size`` values will result in a lot of potentially extra iteration. This parameter can have a strongly negative impact on computation performance. + +* This filter will mark all returns deemed to be ground returns with a + classification value of 2 (per the LAS specification). To extract only these + returns, users can add a :ref:`range filter` to the pipeline. + +.. code-block:: json + + { + "type":"filters.range", + "limits":"Classification[2:2]" + } .. note:: - [Zhang2003]_ describes the consequences and relationships of the - parameters in more detail and is the canonnical resource on the - topic. + + [Zhang2003]_ describes the consequences and relationships of the parameters + in more detail and is the canonnical resource on the topic. Options ------------------------------------------------------------------------------- @@ -75,11 +86,12 @@ initial_distance cell_size Cell Size. [Default: **1**] -classify - Apply classification labels? [Default: **true**] - -extract - Extract ground returns? [Default: **false**] - approximate - Use approximate algorithm? [Default:: **false**] + Use approximate algorithm? [Default: **false**] + +ignore + Optional range of values to ignore. + +last + Consider only last returns (when return information is available)? [Default: + **true**] diff --git a/doc/workshop/exercises/analysis/denoising/denoise.json b/doc/workshop/exercises/analysis/denoising/denoise.json index d9a4d8267b..e166a9967c 100644 --- a/doc/workshop/exercises/analysis/denoising/denoise.json +++ b/doc/workshop/exercises/analysis/denoising/denoise.json @@ -4,13 +4,12 @@ { "type": "filters.outlier", "method": "statistical", - "extract": "true", "multiplier": 3, "mean_k": 8 }, { "type": "filters.range", - "limits": "Z[-100:3000]" + "limits": "Classification![7:7],Z[-100:3000]" }, { "type": "writers.las", diff --git a/doc/workshop/exercises/analysis/denoising/denoising.rst b/doc/workshop/exercises/analysis/denoising/denoising.rst index 8c127500b2..2f70e0bd93 100644 --- a/doc/workshop/exercises/analysis/denoising/denoising.rst +++ b/doc/workshop/exercises/analysis/denoising/denoising.rst @@ -28,7 +28,7 @@ This exercise uses PDAL to remove unwanted noise in an ALS collection. Exercise -------------------------------------------------------------------------------- -PDAL provides a :ref:`filter ` through |PCL| to apply a statistical +PDAL provides the :ref:`outlier filter` to apply a statistical filter to data. Because this operation is somewhat complex, we are going to use a pipeline to @@ -58,14 +58,14 @@ point cloud file we're going to read. 2. :ref:`filters.outlier` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The :ref:`filters.outlier` PDAL filter does most of the work for this operation. +The PDAL :ref:`outlier filter` does most of the work for this +operation. :: { "type": "filters.outlier", "method": "statistical", - "extract": "true", "multiplier": 3, "mean_k": 8 }, @@ -75,17 +75,26 @@ The :ref:`filters.outlier` PDAL filter does most of the work for this operation. 3. :ref:`filters.range` ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +At this point, the outliers have been classified per the LAS specification as +low/noise points with a classification value of 7. The :ref:`range +filter` can remove these noise points by constructing a +:ref:`range ` with the value ``Classification![7:7]``, which passes +every point with a ``Classification`` value **not** equal to 7. + Even with the :ref:`filters.outlier` operation, there is still a cluster of points with extremely negative ``Z`` values. These are some artifact or -miscomputation of processing, and we don't want these points. We are going to -use ::ref:`filters.range` to keep only points that are within the range +miscomputation of processing, and we don't want these points. We can construct +another :ref:`range ` to keep only points that are within the range ``-100 <= Z <= 3000``. +Both :ref:`ranges ` are passed as a comma-separated list to the +:ref:`range filter` via the ``limits`` option. + :: { "type": "filters.range", - "limits": "Z[-100:3000]" + "limits": "Classification![7:7],Z[-100:3000]" }, 4. :ref:`writers.las` diff --git a/doc/workshop/exercises/analysis/ground/ground-run-ground-only.txt b/doc/workshop/exercises/analysis/ground/ground-run-ground-only.txt index e7f6e1be12..cee9c5cc65 100644 --- a/doc/workshop/exercises/analysis/ground/ground-run-ground-only.txt +++ b/doc/workshop/exercises/analysis/ground/ground-run-ground-only.txt @@ -2,5 +2,5 @@ docker run -v /c/Users/Howard/PDAL:/data -t pdal/pdal \ pdal ground \ /data/exercises/analysis/ground/CSite1_orig-utm.laz \ -o /data/exercises/analysis/ground/ground-only.laz \ - --classify=true --extract=true \ + --extract=true \ --writers.las.compression=true --verbose 4 diff --git a/doc/workshop/exercises/analysis/ground/ground.rst b/doc/workshop/exercises/analysis/ground/ground.rst index e18954e5f3..72a5fb5257 100644 --- a/doc/workshop/exercises/analysis/ground/ground.rst +++ b/doc/workshop/exercises/analysis/ground/ground.rst @@ -67,25 +67,28 @@ technique we learned about in :ref:`denoising`. .. literalinclude:: ./ground-run-ground-only.txt :linenos: - :emphasize-lines: 6 + :emphasize-lines: 5 .. note:: - The ``filters.pmf.extract=true`` item causes all data except - ground-classified points to be removed from the set. - -Buildings and other non-ground points are removed with the ``extract`` option -of :ref:`filters.pmf` + The ``--extract=true`` option causes all data except ground-classified + points to be removed from the set. .. image:: ../../../images/ground-ground-only-view.png -2. Now we will remove the noise, using the :ref:`translate_command` to stack the +2. Now we will instead use the :ref:`translate_command` command to stack the :ref:`filters.outlier` and :ref:`filters.pmf` stages: .. literalinclude:: ./translate-run-ground-only.txt :linenos: +In this invocation, we have more control over the process. First the outlier +filter merely classifies outliers with a ``Classification`` value of 7. These +outliers are then ignored during PMF processing with the ``ignore`` option. +Finally, we add a range filter to extract only the ground returns (i.e., +``Classification`` value of 2). + The result is a more accurate representation of the ground returns. .. image:: ../../../images/ground-filtered.png diff --git a/doc/workshop/exercises/analysis/ground/translate-run-ground-only.txt b/doc/workshop/exercises/analysis/ground/translate-run-ground-only.txt index 2602e438c3..b5dc749b51 100644 --- a/doc/workshop/exercises/analysis/ground/translate-run-ground-only.txt +++ b/doc/workshop/exercises/analysis/ground/translate-run-ground-only.txt @@ -2,10 +2,11 @@ docker run -v /c/Users/Howard/PDAL:/data -t pdal/pdal \ pdal translate \ /data/exercises/analysis/ground/CSite1_orig-utm.laz \ -o /data/exercises/analysis/ground/denoised-ground-only.laz \ - outlier pmf \ + outlier pmf range \ --filters.outlier.method="statistical" \ --filters.outlier.mean_k=8 \ --filters.outlier.multiplier=3.0 \ --filters.pmf.cell_size=1.5 \ - --filters.pmf.extract=true \ + --filters.pmf.ignore="Classification[7:7]" \ + --filters.range.limits="Classification[2:2]" \ --writers.las.compression=true --verbose 4 From 284af8a883f146b2652a3bccf029d507de7a22c9 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Mon, 17 Apr 2017 13:32:23 -0500 Subject: [PATCH 11/57] Allow tag names to contain underscores and capital letters. Catch bad tag names during parse. Close #1568 --- pdal/Kernel.cpp | 13 ++++++++----- pdal/PipelineReaderJSON.cpp | 8 +++++--- pdal/Stage.cpp | 28 ++++++++++++++++++++++++++++ pdal/Stage.hpp | 21 +++++++++++++++++++++ test/data/pipeline/options.json.in | 5 +++-- 5 files changed, 65 insertions(+), 10 deletions(-) diff --git a/pdal/Kernel.cpp b/pdal/Kernel.cpp index 1b2ca01cbe..075411a86d 100644 --- a/pdal/Kernel.cpp +++ b/pdal/Kernel.cpp @@ -55,6 +55,8 @@ namespace pdal Kernel::Kernel() : m_showTime(false), m_hardCoreDebug(false) {} + +// Overridden in PipelineKernel to accept "stage" as well. bool Kernel::isStagePrefix(const std::string& stageType) { return (stageType == "readers" || stageType == "writers" || @@ -82,8 +84,6 @@ bool Kernel::parseStageOption(std::string o, std::string& stage, // a bit better than the cast solution. auto islc = [](char c) { return std::islower(c); }; - auto islcOrDigit = [](char c) - { return std::islower(c) || std::isdigit(c); }; std::string::size_type pos = 0; std::string::size_type count = 0; @@ -98,10 +98,13 @@ bool Kernel::parseStageOption(std::string o, std::string& stage, return false; // Get stage_name. - count = Utils::extract(o, pos, islcOrDigit); - if (std::isdigit(o[pos])) + bool ok; + if (stageType == "stage") + ok = Stage::parseTagName(o, pos); + else + ok = Stage::parseName(o, pos); + if (!ok) return false; - pos += count; stage = o.substr(0, pos); if (pos >= o.length() || o[pos++] != '.') return false; diff --git a/pdal/PipelineReaderJSON.cpp b/pdal/PipelineReaderJSON.cpp index 5a5dbb9dd4..d41c74f22c 100644 --- a/pdal/PipelineReaderJSON.cpp +++ b/pdal/PipelineReaderJSON.cpp @@ -252,10 +252,12 @@ std::string PipelineReaderJSON::extractTag(Json::Value& node, TagMap& tags) if (node.isMember("tag")) throw pdal_error("JSON pipeline: found duplicate 'tag' " "entry in stage definition."); + std::string::size_type pos = 0; + if (!Stage::parseTagName(tag, pos) || pos != tag.size()) + throw pdal_error("JSON pipeline: Invalid tag name '" + tag + "'. " + "Must start with letter. Remainder can be letters, " + "digits or underscores."); } - if (Utils::contains(tag, '.')) - throw pdal_error("JSON pipeline: Stage tag name can't contain " - "'.' character."); return tag; } diff --git a/pdal/Stage.cpp b/pdal/Stage.cpp index ac6b1c59c8..76971a483c 100644 --- a/pdal/Stage.cpp +++ b/pdal/Stage.cpp @@ -470,6 +470,34 @@ void Stage::setSpatialReference(MetadataNode& m, } +bool Stage::parseName(std::string o, std::string::size_type& pos) +{ + auto isStageChar = [](char c) + { return std::islower(c) || std::isdigit(c); }; + + std::string::size_type start = pos; + if (!std::islower(o[pos])) + return false; + pos++; + pos += Utils::extract(o, pos, isStageChar); + return true; +} + + +bool Stage::parseTagName(std::string o, std::string::size_type& pos) +{ + auto isTagChar = [](char c) + { return std::isalnum(c) || c == '_'; }; + + std::string::size_type start = pos; + if (!std::isalpha(o[pos])) + return false; + pos++; + pos += Utils::extract(o, pos, isTagChar); + return true; +} + + void Stage::throwError(const std::string& s) const { throw pdal_error(getName() + ": " + s); diff --git a/pdal/Stage.hpp b/pdal/Stage.hpp index 8f0e4473a6..d3d5c3ab03 100644 --- a/pdal/Stage.hpp +++ b/pdal/Stage.hpp @@ -292,6 +292,27 @@ class PDAL_DLL Stage */ void serialize(MetadataNode root, PipelineWriter::TagMap& tags) const; + /** + Parse a stage name from a string. Return the name and update the + position in the input string to the end of the stage name. + + \param o Input string to parse. + \param pos Parsing start/end position. + \return Whether the parsed name is a valid stage name. + */ + static bool parseName(std::string o, std::string::size_type& pos); + + /** + Parse a tag name from a string. Return the name and update the + position in the input string to the end of the tag name. + + \param o Input string to parse. + \param pos Parsing start/end position. + \param tag Parsed tag name. + \return Whether the parsed name is a valid tag name. + */ + static bool parseTagName(std::string o, std::string::size_type& pos); + protected: Options m_options; ///< Stage's options. MetadataNode m_metadata; ///< Stage's metadata. diff --git a/test/data/pipeline/options.json.in b/test/data/pipeline/options.json.in index 90cc3c6cf1..f9a7715b00 100644 --- a/test/data/pipeline/options.json.in +++ b/test/data/pipeline/options.json.in @@ -3,12 +3,13 @@ { "filename" :"@CMAKE_SOURCE_DIR@/test/data/las/1.2-with-color.las", "compression" : "laszip", - "tag": "reader" + "tag": "my_reader" }, { "type":"filters.assign", "assignment":"Z[:]=25", - "tag":"assigner" + "tag":"assigner", + "inputs":"my_reader" }, "@CMAKE_SOURCE_DIR@/test/temp/assigned.las" ] From 2cd676f777101186e5088f0da93c2fdcb96460c2 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Mon, 17 Apr 2017 15:24:32 -0500 Subject: [PATCH 12/57] Hook greedy triangulation into filters. --- filters/GreedyProjection.cpp | 20 ++++++++-- filters/GreedyProjection.hpp | 20 ++++------ pdal/Mesh.hpp | 71 ++++++++++++++++++++++++++++++++++++ pdal/PointView.hpp | 5 +++ pdal/StageFactory.cpp | 2 + 5 files changed, 103 insertions(+), 15 deletions(-) create mode 100644 pdal/Mesh.hpp diff --git a/filters/GreedyProjection.cpp b/filters/GreedyProjection.cpp index b56ca9a540..35ff052659 100644 --- a/filters/GreedyProjection.cpp +++ b/filters/GreedyProjection.cpp @@ -38,12 +38,25 @@ #include #include +#include #include "GreedyProjection.hpp" namespace pdal { +static PluginInfo const s_info = + PluginInfo("filters.greedygrid", "Greedy Triangulation filter", + "http://pdal.io/stages/filters.greedygrid.html"); + +CREATE_STATIC_PLUGIN(1, 0, GreedyProjection, Filter, s_info) + +std::string GreedyProjection::getName() const +{ + return s_info.name; +} + + void GreedyProjection::addArgs(ProgramArgs& args) { args.add("multiplier", "Nearest neighbor distance multiplier", @@ -96,13 +109,14 @@ Eigen::Vector3d GreedyProjection::getNormalCoord(PointId id) void GreedyProjection::addTriangle(PointId a, PointId b, PointId c) { - grid_.emplace_back(a, b, c); + view_->mesh().add(a, b, c); } void GreedyProjection::filter(PointView& view) { KD3Index tree(view); + tree.build(); view_ = &view; const double sqr_mu = mu_ * mu_; @@ -1170,8 +1184,8 @@ void GreedyProjection::filter(PointView& view) } } } - log()->get(LogLevel::Debug) << "Number of triangles: " << grid_.size() << - ".\n"; + log()->get(LogLevel::Debug) << "Number of triangles: " << + view_->mesh().size() << ".\n"; log()->get(LogLevel::Debug) << "Number of unconnected parts: " << nr_parts << ".\n"; if (increase_nnn4fn > 0) diff --git a/filters/GreedyProjection.hpp b/filters/GreedyProjection.hpp index 66e1e59f4f..da615993bb 100644 --- a/filters/GreedyProjection.hpp +++ b/filters/GreedyProjection.hpp @@ -42,8 +42,12 @@ #include #include #include +#include #include +extern "C" int32_t GreedyProjection_ExitFunc(); +extern "C" PF_ExitFunc GreedyProjection_InitPlugin(); + namespace pdal { /** \brief Returns if a point X is visible from point R (or the origin) @@ -177,6 +181,10 @@ namespace pdal view_(nullptr) {}; + static void * create(); + static int32_t destroy(void *); + std::string getName() const; + /** \brief Don't consider points for triangulation if their normal deviates more than this value from the query point's normal. * \param[in] eps_angle maximum surface angle * \note As normal estimation methods usually give smooth transitions at sharp edges, this ensures correct triangulation @@ -270,16 +278,6 @@ namespace pdal Eigen::Vector2d second; }; - struct Triangle - { - Triangle(PointId ta, PointId tb, PointId tc) : a(ta), b(tb), c(tc) - {} - - PointId a; - PointId b; - PointId c; - }; - // Variables made global to decrease the number of parameters to helper functions /** \brief A list of angles to neighbors **/ @@ -338,8 +336,6 @@ namespace pdal Eigen::Vector2d uvn_next_sfn_; /** \brief Temporary variable to store 3 coordiantes **/ Eigen::Vector3d tmp_; - /** \brief Output grid **/ - std::vector grid_; /** \brief Pointer to current point view. **/ PointView *view_; diff --git a/pdal/Mesh.hpp b/pdal/Mesh.hpp new file mode 100644 index 0000000000..eb69def42a --- /dev/null +++ b/pdal/Mesh.hpp @@ -0,0 +1,71 @@ +/****************************************************************************** +* Copyright (c) 2017, Hobu Inc. +* +* 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 + +namespace pdal +{ + +class Triangle +{ +public: + Triangle(PointId a, PointId b, PointId c) : m_a(a), m_b(b), m_c(c) + {} + + PointId m_a; + PointId m_b; + PointId m_c; +}; + +/** + A mesh is a way to represent a set of points connected by edges. Point + indices are into a +*/ +class PDAL_DLL TriangularMesh +{ +public: + TriangularMesh() + {} + + size_t size() const + { return m_index.size(); } + void add(PointId a, PointId b, PointId c) + { m_index.emplace_back(a, b, c); } +protected: + std::deque m_index; +}; + +} // namespace pdal diff --git a/pdal/PointView.hpp b/pdal/PointView.hpp index 3031b0b014..f694372c19 100644 --- a/pdal/PointView.hpp +++ b/pdal/PointView.hpp @@ -36,6 +36,7 @@ #include #include +#include #include #include #include @@ -272,6 +273,9 @@ class PDAL_DLL PointView : public PointContainer } MetadataNode toMetadata() const; + TriangularMesh& mesh() + { return m_mesh; } + protected: PointTableRef m_pointTable; std::deque m_index; @@ -281,6 +285,7 @@ class PDAL_DLL PointView : public PointContainer int m_id; std::queue m_temps; SpatialReference m_spatialReference; + TriangularMesh m_mesh; private: static int m_lastId; diff --git a/pdal/StageFactory.cpp b/pdal/StageFactory.cpp index 676df532ee..f96545f4dc 100644 --- a/pdal/StageFactory.cpp +++ b/pdal/StageFactory.cpp @@ -50,6 +50,7 @@ #include #include #include +#include #include #include #include @@ -251,6 +252,7 @@ StageFactory::StageFactory(bool no_plugins) PluginManager::initializePlugin(ELMFilter_InitPlugin); PluginManager::initializePlugin(EstimateRankFilter_InitPlugin); PluginManager::initializePlugin(FerryFilter_InitPlugin); + PluginManager::initializePlugin(GreedyProjection_InitPlugin); PluginManager::initializePlugin(GroupByFilter_InitPlugin); PluginManager::initializePlugin(HAGFilter_InitPlugin); PluginManager::initializePlugin(IQRFilter_InitPlugin); From b6a4b22b6e6b681e38f83b63a3d0e71080508cab Mon Sep 17 00:00:00 2001 From: Pete Gadomski Date: Mon, 17 Apr 2017 15:52:59 -0600 Subject: [PATCH 13/57] Add a buffer option to the splitter This buffer grows each square by the specified amount, allowing for an overlap between tiles. This is particularily useful for tiled change detection and other such algorithms. --- doc/stages/filters.splitter.rst | 3 ++ filters/SplitterFilter.cpp | 57 +++++++++++++++++++++----- filters/SplitterFilter.hpp | 3 ++ test/unit/filters/SplitterTest.cpp | 64 ++++++++++++++++++++++++++++++ 4 files changed, 116 insertions(+), 11 deletions(-) diff --git a/doc/stages/filters.splitter.rst b/doc/stages/filters.splitter.rst index 6c5fc98456..9943e61c54 100644 --- a/doc/stages/filters.splitter.rst +++ b/doc/stages/filters.splitter.rst @@ -48,3 +48,6 @@ origin_x origin_y Y Origin of the tiles. [Default: none (chosen arbitarily)] +buffer + Amount of overlap to include in each tile. This buffer is added onto length in both the x and the y direction. + [Default: 0.0] diff --git a/filters/SplitterFilter.cpp b/filters/SplitterFilter.cpp index 52608c6136..e390ecd9e1 100644 --- a/filters/SplitterFilter.cpp +++ b/filters/SplitterFilter.cpp @@ -63,8 +63,18 @@ void SplitterFilter::addArgs(ProgramArgs& args) std::numeric_limits::quiet_NaN()); args.add("origin_y", "Y origin for a cell", m_yOrigin, std::numeric_limits::quiet_NaN()); + args.add("buffer", "Size of buffer (overlap) to include around each tile.", + m_buffer, 0.0); } +void SplitterFilter::initialize() { + if (!(m_buffer < m_length / 2.)) { + std::stringstream oss; + oss << "Buffer (" << m_buffer << + ") must be less than half of length (" << m_length << ")"; + throw pdal_error(oss.str()); + } +} PointViewSet SplitterFilter::run(PointViewPtr inView) { @@ -72,6 +82,14 @@ PointViewSet SplitterFilter::run(PointViewPtr inView) if (!inView->size()) return viewSet; + auto addPoint = [this, &inView](PointId idx, int xpos, int ypos) { + Coord loc(xpos, ypos); + PointViewPtr& outView = m_viewMap[loc]; + if (!outView) + outView = inView->makeNew(); + outView->appendPoint(*inView.get(), idx); + }; + // Use the location of the first point as the origin, unless specified. // (!= test == isnan(), which doesn't exist on windows) if (m_xOrigin != m_xOrigin) @@ -84,22 +102,31 @@ PointViewSet SplitterFilter::run(PointViewPtr inView) for (PointId idx = 0; idx < inView->size(); idx++) { double x = inView->getFieldAs(Dimension::Id::X, idx); - x -= m_xOrigin; - int xpos = x / m_length; - if (x < 0) + double dx = x - m_xOrigin; + int xpos = dx / m_length; + if (dx < 0) xpos--; double y = inView->getFieldAs(Dimension::Id::Y, idx); - y -= m_yOrigin; - int ypos = y / m_length; - if (y < 0) + double dy = y - m_yOrigin; + int ypos = dy / m_length; + if (dy < 0) ypos--; - Coord loc(xpos, ypos); - PointViewPtr& outView = m_viewMap[loc]; - if (!outView) - outView = inView->makeNew(); - outView->appendPoint(*inView.get(), idx); + addPoint(idx, xpos, ypos); + + if (m_buffer > 0.0) { + if (squareContains(xpos - 1, ypos, x, y)) { + addPoint(idx, xpos - 1, ypos); + } else if (squareContains(xpos + 1, ypos, x, y)) { + addPoint(idx, xpos + 1, ypos); + } + if (squareContains(xpos, ypos - 1, x, y)) { + addPoint(idx, xpos, ypos - 1); + } else if (squareContains(xpos, ypos + 1, x, y)) { + addPoint(idx, xpos, ypos + 1); + } + } } // Pull the buffers out of the map and stick them in the standard @@ -109,4 +136,12 @@ PointViewSet SplitterFilter::run(PointViewPtr inView) return viewSet; } +bool SplitterFilter::squareContains(int xpos, int ypos, double x, double y) const { + double minx = m_xOrigin + xpos * m_length - m_buffer; + double maxx = minx + m_length + 2 * m_buffer; + double miny = m_yOrigin + ypos * m_length - m_buffer; + double maxy = miny + m_length + 2 * m_buffer; + return minx < x && x < maxx && miny < y && y < maxy; +} + } // pdal diff --git a/filters/SplitterFilter.hpp b/filters/SplitterFilter.hpp index 6711b35016..8c79197f20 100644 --- a/filters/SplitterFilter.hpp +++ b/filters/SplitterFilter.hpp @@ -71,10 +71,13 @@ class PDAL_DLL SplitterFilter : public pdal::Filter double m_length; double m_xOrigin; double m_yOrigin; + double m_buffer; std::map m_viewMap; virtual void addArgs(ProgramArgs& args); + virtual void initialize(); virtual PointViewSet run(PointViewPtr view); + bool squareContains(int xpos, int ypos, double x, double y) const; SplitterFilter& operator=(const SplitterFilter&); // not implemented SplitterFilter(const SplitterFilter&); // not implemented diff --git a/test/unit/filters/SplitterTest.cpp b/test/unit/filters/SplitterTest.cpp index 6177093e58..9bc7c5e9ea 100644 --- a/test/unit/filters/SplitterTest.cpp +++ b/test/unit/filters/SplitterTest.cpp @@ -113,3 +113,67 @@ TEST(SplitterTest, test_tile_filter) EXPECT_EQ(view->size(), counts[i]); } } + +TEST(SplitterTest, test_buffer) +{ + Options readerOptions; + readerOptions.add("filename", Support::datapath("las/1.2-with-color.las")); + LasReader reader; + reader.setOptions(readerOptions); + + Options splitterOptions; + splitterOptions.add("length", 1000); + splitterOptions.add("buffer", 20); + + SplitterFilter splitter; + splitter.setOptions(splitterOptions); + splitter.setInput(reader); + + PointTable table; + PointViewPtr view(new PointView(table)); + splitter.prepare(table); + + StageWrapper::ready(reader, table); + PointViewSet viewSet = StageWrapper::run(reader, view); + StageWrapper::done(reader, table); + ASSERT_EQ(viewSet.size(), 1u); + view = *viewSet.begin(); + + StageWrapper::ready(splitter, table); + viewSet = StageWrapper::run(splitter, view); + StageWrapper::done(splitter, table); + + std::vector views; + std::map bounds; + + for (auto it = viewSet.begin(); it != viewSet.end(); ++it) + { + BOX2D b; + PointViewPtr v = *it; + v->calculateBounds(b); + EXPECT_TRUE(b.maxx - b.minx <= 1040); + EXPECT_TRUE(b.maxy - b.miny <= 1040); + bounds[v] = b; + views.push_back(v); + } + + auto sorter = [&bounds](PointViewPtr p1, PointViewPtr p2) + { + BOX2D b1 = bounds[p1]; + BOX2D b2 = bounds[p2]; + + return b1.minx < b2.minx ? true : + b1.minx > b2.minx ? false : + b1.miny < b2.miny; + }; + std::sort(views.begin(), views.end(), sorter); + + EXPECT_EQ(views.size(), 24u); + size_t counts[] = {26, 26, 3, 28, 27, 13, 65, 80, 47, 80, 89, 12, 94, + 77, 5, 79, 65, 34, 63, 67, 74, 69, 36, 5}; + for (size_t i = 0; i < views.size(); ++i) + { + PointViewPtr view = views[i]; + EXPECT_EQ(view->size(), counts[i]); + } +} From 82a09ecdcd59db917c872dfa7a27ddaa7504ad8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Lemoine?= Date: Tue, 18 Apr 2017 10:06:36 +0200 Subject: [PATCH 14/57] Add note to docs about pgpointcloud tests --- doc/development/testing.rst | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/doc/development/testing.rst b/doc/development/testing.rst index d3cab6aee0..7997a1cffb 100644 --- a/doc/development/testing.rst +++ b/doc/development/testing.rst @@ -312,6 +312,18 @@ the available flags type:: Key among these flags are the ability to list tests (``--gtest_list_tests``) and to run only select tests (``--gtest_filter``). +.. note:: + + If the PostgreSQL PointCloud plugin was enabled on the CMake command line + (with ``-DBUILD_PLUGIN_PGPOINTCLOUD=ON``) then `ctest` will attempt to run + the ``pgpointcloud`` tests. And you will get PostgreSQL connection errors + if the `libpq environment variables`_ are not correctly set in your shell. + This is for example how you can run the ``pgpointcloud`` tests:: + + $ PGUSER=pdal PGPASSWORD=pdal PGHOST=localhost ctest -R pgpointcloudtest + +.. _`libpg environment variables`: https://www.postgresql.org/docs/current/static/libpq-envars.html + Test Data ========= From 99f8359cd892b220815cfdc7d68ac380b1bd6e57 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Tue, 18 Apr 2017 11:51:18 -0500 Subject: [PATCH 15/57] Hooks for multiple meshes. --- filters/GreedyProjection.cpp | 5 +++-- filters/GreedyProjection.hpp | 7 ++++++- pdal/Mesh.hpp | 11 +++++++++-- pdal/PointView.cpp | 27 +++++++++++++++++++++++++++ pdal/PointView.hpp | 23 ++++++++++++++++++----- 5 files changed, 63 insertions(+), 10 deletions(-) diff --git a/filters/GreedyProjection.cpp b/filters/GreedyProjection.cpp index 35ff052659..a0df302a05 100644 --- a/filters/GreedyProjection.cpp +++ b/filters/GreedyProjection.cpp @@ -109,7 +109,7 @@ Eigen::Vector3d GreedyProjection::getNormalCoord(PointId id) void GreedyProjection::addTriangle(PointId a, PointId b, PointId c) { - view_->mesh().add(a, b, c); + mesh_->add(a, b, c); } @@ -119,6 +119,7 @@ void GreedyProjection::filter(PointView& view) tree.build(); view_ = &view; + mesh_ = view_->createMesh(getName()); const double sqr_mu = mu_ * mu_; const double sqr_max_edge = search_radius_*search_radius_; @@ -1185,7 +1186,7 @@ void GreedyProjection::filter(PointView& view) } } log()->get(LogLevel::Debug) << "Number of triangles: " << - view_->mesh().size() << ".\n"; + mesh_->size() << ".\n"; log()->get(LogLevel::Debug) << "Number of unconnected parts: " << nr_parts << ".\n"; if (increase_nnn4fn > 0) diff --git a/filters/GreedyProjection.hpp b/filters/GreedyProjection.hpp index da615993bb..d722a86f10 100644 --- a/filters/GreedyProjection.hpp +++ b/filters/GreedyProjection.hpp @@ -50,6 +50,8 @@ extern "C" PF_ExitFunc GreedyProjection_InitPlugin(); namespace pdal { + class TriangularMesh; + /** \brief Returns if a point X is visible from point R (or the origin) * when taking into account the segment between the points S1 and S2 * \param X 2D coordinate of the point @@ -178,7 +180,8 @@ namespace pdal uvn_next_ffn_ (), uvn_next_sfn_ (), tmp_ (), - view_(nullptr) + view_(nullptr), + mesh_(nullptr) {}; static void * create(); @@ -338,6 +341,8 @@ namespace pdal Eigen::Vector3d tmp_; /** \brief Pointer to current point view. **/ PointView *view_; + /** \brief Pointer to the mesh we're creating. **/ + TriangularMesh *mesh_; /** \brief Forms a new triangle by connecting the current neighbor to the query point * and the previous neighbor diff --git a/pdal/Mesh.hpp b/pdal/Mesh.hpp index eb69def42a..afe26d228d 100644 --- a/pdal/Mesh.hpp +++ b/pdal/Mesh.hpp @@ -52,9 +52,16 @@ class Triangle /** A mesh is a way to represent a set of points connected by edges. Point - indices are into a + indices are into a point view. */ -class PDAL_DLL TriangularMesh +class PDAL_DLL Mesh +{}; + + +/** + A mesh where the faces are triangles. +*/ +class PDAL_DLL TriangularMesh : public Mesh { public: TriangularMesh() diff --git a/pdal/PointView.cpp b/pdal/PointView.cpp index 36b542e5aa..d3574b3ed7 100644 --- a/pdal/PointView.cpp +++ b/pdal/PointView.cpp @@ -36,6 +36,7 @@ #include #include +#include namespace pdal { @@ -54,6 +55,12 @@ PointView::PointView(PointTableRef pointTable, const SpatialReference& srs) : m_id = ++m_lastId; } +PointView::~PointView() +{ + for (auto it = m_meshes.begin(); it != m_meshes.end(); ++it) + delete it->second; +} + PointViewIter PointView::begin() { return PointViewIter(this, 0); @@ -155,6 +162,26 @@ MetadataNode PointView::toMetadata() const } +TriangularMesh *PointView::createMesh(const std::string& name) +{ + if (Utils::contains(m_meshes, name)) + return nullptr; + auto res = m_meshes.insert(std::make_pair(name, new TriangularMesh)); + if (res.second) + return res.first->second; + return nullptr; +} + + +TriangularMesh *PointView::mesh(const std::string& name) +{ + auto it = m_meshes.find(name); + if (it != m_meshes.end()) + return it->second; + return nullptr; +} + + void PointView::dump(std::ostream& ostr) const { using std::endl; diff --git a/pdal/PointView.hpp b/pdal/PointView.hpp index f694372c19..4fd4ed0752 100644 --- a/pdal/PointView.hpp +++ b/pdal/PointView.hpp @@ -76,8 +76,7 @@ class PDAL_DLL PointView : public PointContainer PointView(PointTableRef pointTable); PointView(PointTableRef pointTable, const SpatialReference& srs); - virtual ~PointView() - {} + virtual ~PointView(); PointViewIter begin(); PointViewIter end(); @@ -273,8 +272,22 @@ class PDAL_DLL PointView : public PointContainer } MetadataNode toMetadata() const; - TriangularMesh& mesh() - { return m_mesh; } + /** + Creates a mesh with the specified name. + + \param name Name of the mesh. + \return Pointer to the new mesh. Null is returned if the mesh + already exists. + */ + TriangularMesh *createMesh(const std::string& name); + + /** + Get a pointer to a mesh. + + \param name Name of the mesh. + \return New mesh. Null is returned if the mesh already exists. + */ + TriangularMesh *mesh(const std::string& name); protected: PointTableRef m_pointTable; @@ -285,7 +298,7 @@ class PDAL_DLL PointView : public PointContainer int m_id; std::queue m_temps; SpatialReference m_spatialReference; - TriangularMesh m_mesh; + std::map m_meshes; private: static int m_lastId; From 05b79909eb00a55f0048d951f396be2275140352 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Wed, 19 Apr 2017 13:11:42 -0500 Subject: [PATCH 16/57] Remove MultiFilter. --- filters/MergeFilter.hpp | 2 +- pdal/Filter.hpp | 7 ------- 2 files changed, 1 insertion(+), 8 deletions(-) diff --git a/filters/MergeFilter.hpp b/filters/MergeFilter.hpp index a73b2c221e..f8a2f3fb34 100644 --- a/filters/MergeFilter.hpp +++ b/filters/MergeFilter.hpp @@ -43,7 +43,7 @@ extern "C" PF_ExitFunc MergeFilter_InitPlugin(); namespace pdal { -class PDAL_DLL MergeFilter : public MultiFilter +class PDAL_DLL MergeFilter : public Filter { public: MergeFilter () diff --git a/pdal/Filter.hpp b/pdal/Filter.hpp index 7ac159d2ad..35be2b65db 100644 --- a/pdal/Filter.hpp +++ b/pdal/Filter.hpp @@ -65,12 +65,5 @@ class PDAL_DLL Filter : public Stage Filter(const Filter&); // not implemented }; -class PDAL_DLL MultiFilter : public Filter -{ -public: - MultiFilter() : Filter() - {} -}; - } // namespace pdal From e9ad7088155aeb26a7b5e6e167aa93d5bd061ce6 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Wed, 19 Apr 2017 15:04:30 -0500 Subject: [PATCH 17/57] Enhance PLY writer. Close #1572 --- io/PlyWriter.cpp | 72 ++++++++++++++++++++--- io/PlyWriter.hpp | 6 ++ pdal/Mesh.hpp | 80 +++++++++++++++++++++++++ pdal/PointView.cpp | 30 ++++++++++ pdal/PointView.hpp | 26 ++++++-- test/unit/PointViewTest.cpp | 114 ------------------------------------ 6 files changed, 203 insertions(+), 125 deletions(-) create mode 100644 pdal/Mesh.hpp diff --git a/io/PlyWriter.cpp b/io/PlyWriter.cpp index 407a8e1a05..89564276d8 100644 --- a/io/PlyWriter.cpp +++ b/io/PlyWriter.cpp @@ -34,6 +34,7 @@ #include "PlyWriter.hpp" +#include #include #include @@ -92,6 +93,8 @@ void PlyWriter::addArgs(ProgramArgs& args) { args.add("filename", "Output filename", m_filename).setPositional(); args.add("storage_mode", "PLY Storage mode", m_storageModeSpec, "default"); + args.add("dims", "Dimension names", m_dimNames); + args.add("faces", "Write faces", m_faces); } @@ -125,6 +128,28 @@ void PlyWriter::initialize() } +void PlyWriter::prepared(PointTableRef table) +{ + if (m_dimNames.size()) + { + for (auto& name : m_dimNames) + { + auto id = table.layout()->findDim(name); + if (id == Dimension::Id::Unknown) + throwError("Unknown dimension '" + name + "' in provided " + "dimension list."); + m_dims.push_back(id); + } + } + else + { + m_dims = table.layout()->dims(); + for (auto dim : m_dims) + m_dimNames.push_back(Utils::tolower(table.layout()->dimName(dim))); + } +} + + void PlyWriter::ready(PointTableRef table) { try @@ -144,7 +169,20 @@ void PlyWriter::ready(PointTableRef table) void PlyWriter::write(const PointViewPtr data) { + if (m_faces && m_pointCollector->size()) + throwError("Can't output faces for pipeline with multiple " + "point views."); m_pointCollector->append(*data); + if (m_faces) + { + m_mesh = data->mesh(); + if (!m_mesh) + throwError("Can't write mesh faces. No mesh has been generated."); + } + if (data->size() > std::numeric_limits::max()) + throwError("Can't write PLY file. Only " + + std::to_string(std::numeric_limits::max()) + + " points supported."); } @@ -155,13 +193,19 @@ void PlyWriter::done(PointTableRef table) if (!ply_add_element(m_ply, "vertex", m_pointCollector->size())) throwError("Could not add vertex element"); - auto dimensions = table.layout()->dims(); - for (auto dim : dimensions) { - std::string name = table.layout()->dimName(dim); + auto ni = m_dimNames.begin(); + for (auto dim : m_dims) { + std::string name = *ni++; e_ply_type plyType = getPlyType(table.layout()->dimType(dim)); if (!ply_add_scalar_property(m_ply, name.c_str(), plyType)) throwError("Could not add scalar property '" + name + "'"); } + if (m_faces) + { + ply_add_element(m_ply, "face", m_mesh->size()); + ply_add_list_property(m_ply, "vertex_indices", PLY_UINT8, PLY_UINT32); + } + if (!ply_add_comment(m_ply, "Generated by PDAL")) throwError("Could not add comment"); @@ -170,12 +214,26 @@ void PlyWriter::done(PointTableRef table) for (PointId index = 0; index < m_pointCollector->size(); ++index) { - for (auto dim : dimensions) + auto ni = m_dimNames.begin(); + auto di = m_dims.begin(); + for (; di != m_dims.end(); ++di, ++ni) { - double value = m_pointCollector->getFieldAs(dim, index); + double value = m_pointCollector->getFieldAs(*di, index); if (!ply_write(m_ply, value)) - throwError("Error writing dimension '" + - table.layout()->dimName(dim) + "' of point number " + + throwError("Error writing dimension '" + *ni + "' of point " + "number " + Utils::toString(index) + "."); + } + } + if (m_faces) + { + for (PointId id = 0; id < m_mesh->size(); id++) + { + const Triangle& t = (*m_mesh)[id]; + if (!(ply_write(m_ply, 3) && + ply_write(m_ply, t.m_a) && + ply_write(m_ply, t.m_b) && + ply_write(m_ply, t.m_c))) + throwError("Error writing face number " + Utils::toString(index) + "."); } } diff --git a/io/PlyWriter.hpp b/io/PlyWriter.hpp index 39dff77593..f8b8fc6a30 100644 --- a/io/PlyWriter.hpp +++ b/io/PlyWriter.hpp @@ -43,6 +43,7 @@ extern "C" PF_ExitFunc PlyWriter_InitPlugin(); namespace pdal { +class TriangularMesh; class PDAL_DLL PlyWriter : public Writer { @@ -63,6 +64,7 @@ class PDAL_DLL PlyWriter : public Writer private: virtual void addArgs(ProgramArgs& args); virtual void initialize(); + virtual void prepared(PointTableRef table); virtual void ready(PointTableRef table); virtual void write(const PointViewPtr data); virtual void done(PointTableRef table); @@ -72,6 +74,10 @@ class PDAL_DLL PlyWriter : public Writer PointViewPtr m_pointCollector; std::string m_storageModeSpec; e_ply_storage_mode m_storageMode; + bool m_faces; + StringList m_dimNames; + Dimension::IdList m_dims; + TriangularMesh *m_mesh; }; } diff --git a/pdal/Mesh.hpp b/pdal/Mesh.hpp new file mode 100644 index 0000000000..41e1db4af0 --- /dev/null +++ b/pdal/Mesh.hpp @@ -0,0 +1,80 @@ +/****************************************************************************** +* Copyright (c) 2017, Hobu Inc. +* +* 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 + +namespace pdal +{ + +class Triangle +{ +public: + Triangle(PointId a, PointId b, PointId c) : m_a(a), m_b(b), m_c(c) + {} + + PointId m_a; + PointId m_b; + PointId m_c; +}; + +/** + A mesh is a way to represent a set of points connected by edges. Point + indices are into a point view. +*/ +class PDAL_DLL Mesh +{}; + + +/** + A mesh where the faces are triangles. +*/ +class PDAL_DLL TriangularMesh : public Mesh +{ +public: + TriangularMesh() + {} + + size_t size() const + { return m_index.size(); } + void add(PointId a, PointId b, PointId c) + { m_index.emplace_back(a, b, c); } + const Triangle& operator[](PointId id) const + { return m_index[id]; } +protected: + std::deque m_index; +}; + +} // namespace pdal diff --git a/pdal/PointView.cpp b/pdal/PointView.cpp index 36b542e5aa..3cd446bf83 100644 --- a/pdal/PointView.cpp +++ b/pdal/PointView.cpp @@ -34,8 +34,10 @@ #include +#include #include #include +#include namespace pdal { @@ -54,6 +56,11 @@ PointView::PointView(PointTableRef pointTable, const SpatialReference& srs) : m_id = ++m_lastId; } + +PointView::~PointView() +{} + + PointViewIter PointView::begin() { return PointViewIter(this, 0); @@ -155,6 +162,29 @@ MetadataNode PointView::toMetadata() const } +TriangularMesh *PointView::createMesh(const std::string& name) +{ + if (Utils::contains(m_meshes, name)) + return nullptr; + auto res = m_meshes.insert(std::make_pair(name, + std::unique_ptr(new TriangularMesh))); + if (res.second) + return res.first->second.get(); + return nullptr; +} + + +TriangularMesh *PointView::mesh(const std::string& name) +{ + auto it = m_meshes.find(name); + if (it != m_meshes.end()) + return it->second.get(); + if (name.empty() && m_meshes.size()) + return m_meshes.begin()->second.get(); + return nullptr; +} + + void PointView::dump(std::ostream& ostr) const { using std::endl; diff --git a/pdal/PointView.hpp b/pdal/PointView.hpp index 75687c7696..0aa38353d4 100644 --- a/pdal/PointView.hpp +++ b/pdal/PointView.hpp @@ -36,6 +36,7 @@ #include #include +#include #include #include #include @@ -72,11 +73,11 @@ class PDAL_DLL PointView : public PointContainer friend class PointIdxRef; friend struct PointViewLess; public: + PointView(const PointView&) = delete; + PointView& operator=(const PointView&) = delete; PointView(PointTableRef pointTable); PointView(PointTableRef pointTable, const SpatialReference& srs); - - virtual ~PointView() - {} + virtual ~PointView(); PointViewIter begin(); PointViewIter end(); @@ -243,7 +244,6 @@ class PDAL_DLL PointView : public PointContainer } } - /// Provides access to the memory storing the point data. Though this /// function is public, other access methods are safer and preferred. char *getPoint(PointId id) @@ -273,6 +273,23 @@ class PDAL_DLL PointView : public PointContainer } MetadataNode toMetadata() const; + /** + Creates a mesh with the specified name. + + \param name Name of the mesh. + \return Pointer to the new mesh. Null is returned if the mesh + already exists. + */ + TriangularMesh *createMesh(const std::string& name); + + /** + Get a pointer to a mesh. + + \param name Name of the mesh. + \return New mesh. Null is returned if the mesh already exists. + */ + TriangularMesh *mesh(const std::string& name = ""); + protected: PointTableRef m_pointTable; std::deque m_index; @@ -282,6 +299,7 @@ class PDAL_DLL PointView : public PointContainer int m_id; std::queue m_temps; SpatialReference m_spatialReference; + std::map> m_meshes; private: static int m_lastId; diff --git a/test/unit/PointViewTest.cpp b/test/unit/PointViewTest.cpp index 16b9359ebf..f9b9eb8504 100644 --- a/test/unit/PointViewTest.cpp +++ b/test/unit/PointViewTest.cpp @@ -161,58 +161,6 @@ TEST(PointViewTest, getFloat) } -TEST(PointViewTest, copy) -{ - PointTable table; - PointViewPtr view = makeTestView(table); - - PointView d2(*view); - - // read the view back out - { - EXPECT_EQ( - d2.getFieldAs(Dimension::Id::Classification, 0), - view->getFieldAs(Dimension::Id::Classification, 0)); - EXPECT_EQ(d2.getFieldAs(Dimension::Id::X, 0), - view->getFieldAs(Dimension::Id::X, 0)); - EXPECT_FLOAT_EQ(d2.getFieldAs(Dimension::Id::Y, 0), - view->getFieldAs(Dimension::Id::Y, 0)); - } - - for (int i = 1; i < 17; i++) - { - uint8_t x = d2.getFieldAs(Dimension::Id::Classification, i); - int32_t y = d2.getFieldAs(Dimension::Id::X, i); - double z = d2.getFieldAs(Dimension::Id::Y, i); - - EXPECT_EQ(x, i + 1u); - EXPECT_EQ(y, i * 10); - EXPECT_TRUE(Utils::compare_approx(z, i * 100.0, - (std::numeric_limits::min)())); - } - EXPECT_EQ(view->size(), d2.size()); - -} - -TEST(PointViewTest, copyCtor) -{ - PointTable table; - PointViewPtr view = makeTestView(table); - - PointView d2(*view); - verifyTestView(d2); -} - -TEST(PointViewTest, assignment) -{ - PointTable table; - PointViewPtr view = makeTestView(table); - - PointView d2 = *view; - verifyTestView(d2); -} - - TEST(PointViewTest, bigfile) { PointTable table; @@ -280,68 +228,6 @@ TEST(PointViewTest, bigfile) } -//ABELL - Move to KdIndex -/** -TEST(PointViewTest, kdindex) -{ - LasReader reader(Support::viewpath("1.2-with-color.las")); - reader.prepare(); - - const Schema& schema = reader.getSchema(); - uint32_t capacity(1000); - PointView view(schema, capacity); - - StageSequentialIterator* iter = reader.createSequentialIterator(view); - - { - uint32_t numRead = iter->read(view); - EXPECT_EQ(numRead, capacity); - } - - EXPECT_EQ(view.getCapacity(), capacity); - EXPECT_EQ(view.getSchema(), schema); - - - IndexedPointView iview(view); - EXPECT_EQ(iview.getCapacity(), capacity); - EXPECT_EQ(iview.getSchema(), schema); - - iview.build(); - - unsigned k = 8; - - // If the query distance is 0, just return the k nearest neighbors - std::vector ids = idata.neighbors(636199, 849238, 428.05, k); - EXPECT_EQ(ids.size(), k); - EXPECT_EQ(ids[0], 8u); - EXPECT_EQ(ids[1], 7u); - EXPECT_EQ(ids[2], 9u); - EXPECT_EQ(ids[3], 42u); - EXPECT_EQ(ids[4], 40u); - - std::vector dist_ids = idata.neighbors(636199, 849238, 428.05, 3); - - EXPECT_EQ(dist_ids.size(), 3u); - EXPECT_EQ(dist_ids[0], 8u); - - std::vector nids = idata.neighbors(636199, 849238, 428.05, k); - - EXPECT_EQ(nids.size(), k); - EXPECT_EQ(nids[0], 8u); - EXPECT_EQ(nids[1], 7u); - EXPECT_EQ(nids[2], 9u); - EXPECT_EQ(nids[3], 42u); - EXPECT_EQ(nids[4], 40u); - - std::vector rids = iview.radius(637012.24, 849028.31, - 431.66, 100000); - EXPECT_EQ(rids.size(), 11u); - - delete iter; -} -**/ - - static void check_bounds(const BOX3D& box, double minx, double maxx, double miny, double maxy, From 80284b30c1409120954b1fdadb4ee264ad1dd4dc Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Wed, 19 Apr 2017 15:29:36 -0500 Subject: [PATCH 18/57] Correct index -> id. --- io/PlyWriter.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/io/PlyWriter.cpp b/io/PlyWriter.cpp index 89564276d8..278aae28ee 100644 --- a/io/PlyWriter.cpp +++ b/io/PlyWriter.cpp @@ -234,7 +234,7 @@ void PlyWriter::done(PointTableRef table) ply_write(m_ply, t.m_b) && ply_write(m_ply, t.m_c))) throwError("Error writing face number " + - Utils::toString(index) + "."); + Utils::toString(id) + "."); } } From b1fed0114a81596fd49ae6f734ec04183a3fe84a Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Wed, 19 Apr 2017 16:42:35 -0500 Subject: [PATCH 19/57] Working greedy projection filter. --- filters/GreedyProjection.cpp | 32 ++++++---- filters/GreedyProjection.hpp | 1 + filters/NormalFilter.cpp | 33 +++++----- filters/NormalFilter.hpp | 5 +- io/PlyWriter.cpp | 72 +++++++++++++++++++--- io/PlyWriter.hpp | 6 ++ pdal/Dimension.json | 5 ++ pdal/Mesh.hpp | 2 + pdal/PointView.cpp | 44 +++++++++++--- pdal/PointView.hpp | 14 ++++- test/unit/PointViewTest.cpp | 114 ----------------------------------- 11 files changed, 167 insertions(+), 161 deletions(-) diff --git a/filters/GreedyProjection.cpp b/filters/GreedyProjection.cpp index a0df302a05..dc5b051981 100644 --- a/filters/GreedyProjection.cpp +++ b/filters/GreedyProjection.cpp @@ -39,6 +39,7 @@ #include #include +#include #include "GreedyProjection.hpp" @@ -46,8 +47,8 @@ namespace pdal { static PluginInfo const s_info = - PluginInfo("filters.greedygrid", "Greedy Triangulation filter", - "http://pdal.io/stages/filters.greedygrid.html"); + PluginInfo("filters.greedymesh", "Greedy Triangulation filter", + "http://pdal.io/stages/filters.greedymesh.html"); CREATE_STATIC_PLUGIN(1, 0, GreedyProjection, Filter, s_info) @@ -73,6 +74,14 @@ void GreedyProjection::addArgs(ProgramArgs& args) "consideration", eps_angle_, M_PI / 4); } + +void GreedyProjection::addDimensions(PointLayoutPtr layout) +{ + layout->registerDims( { Dimension::Id::NormalX, Dimension::Id::NormalY, + Dimension::Id::NormalZ } ); +} + + void GreedyProjection::initialize() { if (search_radius_ <= 0) @@ -115,8 +124,9 @@ void GreedyProjection::addTriangle(PointId a, PointId b, PointId c) void GreedyProjection::filter(PointView& view) { - KD3Index tree(view); - tree.build(); + NormalFilter().doFilter(view); + + KD3Index& tree = view.build3dIndex(); view_ = &view; mesh_ = view_->createMesh(getName()); @@ -184,7 +194,7 @@ void GreedyProjection::filter(PointView& view) // Initializing PointId isFree = 0; bool done = false; - int is_free=0, nr_parts=0, increase_nnn4fn=0, increase_nnn4s=0, increase_dist=0, nr_touched = 0; + int nr_parts=0, increase_nnn4fn=0, increase_nnn4s=0, increase_dist=0, nr_touched = 0; bool is_fringe; angles_.resize(nnn_); std::vector > uvn_nn (nnn_); @@ -193,7 +203,7 @@ void GreedyProjection::filter(PointView& view) // iterating through fringe points and finishing them until everything is done while (!done) { - R_ = is_free; + R_ = isFree; if (state_[R_] == GP3Type::FREE) { state_[R_] = GP3Type::NONE; @@ -202,6 +212,7 @@ void GreedyProjection::filter(PointView& view) // creating starting triangle tree.knnSearch(R_, nnn_, &nnIdx, &sqrDists); + double sqr_dist_threshold = (std::min)(sqr_max_edge, sqr_mu * sqrDists[1]); @@ -232,7 +243,6 @@ void GreedyProjection::filter(PointView& view) size_t nr_edge = 0; std::vector doubleEdges; // nearest neighbor with index 0 is the query point R_ itself - //ABELL - Need to check that the above is true. for (int i = 1; i < nnn_; i++) { // Transforming coordinates @@ -271,7 +281,6 @@ void GreedyProjection::filter(PointView& view) // Verify the visibility of each potential new vertex // nearest neighbor with index 0 is the query point R_ itself - // ABELL - Check the above. for (int i = 1; i < nnn_; i++) if ((angles_[i].visible) && (ffn_[R_] != nnIdx[i]) && (sfn_[R_] != nnIdx[i])) @@ -342,7 +351,7 @@ void GreedyProjection::filter(PointView& view) while (not_found); } - // Set the index of the first free point in is_free. + // Set the index of the first free point in isFree. done = true; auto it = std::find(state_.begin(), state_.end(), GP3Type::FREE); if (it == state_.end()) @@ -350,7 +359,7 @@ void GreedyProjection::filter(PointView& view) else { done = false; - is_free = std::distance(state_.begin(), it); + isFree = std::distance(state_.begin(), it); } is_fringe = true; @@ -433,7 +442,6 @@ void GreedyProjection::filter(PointView& view) size_t nr_edge = 0; std::vector doubleEdges; // nearest neighbor with index 0 is the query point R_ itself - //ABELL - Check this. for (int i = 1; i < nnn_; i++) { tmp_ = getCoord(nnIdx[i]) - proj_qp_; @@ -543,7 +551,6 @@ void GreedyProjection::filter(PointView& view) // Verify the visibility of each potential new vertex // nearest neighbor with index 0 is the query point R_ itself - //ABELL - Check above (AGAIN?!) for (int i = 1; i < nnn_; i++) if ((angles_[i].visible) && (ffn_[R_] != nnIdx[i]) && (sfn_[R_] != nnIdx[i])) @@ -714,7 +721,6 @@ void GreedyProjection::filter(PointView& view) int nnCB(-1); // nearest neighbor with index 0 is the query point R_ itself - // ABELL - Again. for (int i = 1; i < nnn_; i++) { // NOTE: nnCB is an index in nnIdx diff --git a/filters/GreedyProjection.hpp b/filters/GreedyProjection.hpp index d722a86f10..61e60a4351 100644 --- a/filters/GreedyProjection.hpp +++ b/filters/GreedyProjection.hpp @@ -433,6 +433,7 @@ namespace pdal } void addArgs(ProgramArgs& args); + void addDimensions(PointLayoutPtr layout); void initialize(); void filter(PointView& view); void addTriangle(PointId a, PointId b, PointId c); diff --git a/filters/NormalFilter.cpp b/filters/NormalFilter.cpp index 4cc4770d3f..35472a42a6 100644 --- a/filters/NormalFilter.cpp +++ b/filters/NormalFilter.cpp @@ -67,19 +67,21 @@ void NormalFilter::addArgs(ProgramArgs& args) void NormalFilter::addDimensions(PointLayoutPtr layout) { - m_nx = layout->registerOrAssignDim("NormalX", Dimension::Type::Double); - m_ny = layout->registerOrAssignDim("NormalY", Dimension::Type::Double); - m_nz = layout->registerOrAssignDim("NormalZ", Dimension::Type::Double); - m_curvature = layout->registerOrAssignDim("Curvature", - Dimension::Type::Double); + layout->registerDims( { Dimension::Id::NormalX, Dimension::Id::NormalY, + Dimension::Id::NormalZ, Dimension::Id::Curvature } ); } -void NormalFilter::filter(PointView& view) + +void NormalFilter::doFilter(PointView& view) { - using namespace Eigen; + m_knn = 8; + filter(view); +} - KD3Index kdi(view); - kdi.build(); + +void NormalFilter::filter(PointView& view) +{ + KD3Index& kdi = view.build3dIndex(); for (PointId i = 0; i < view.size(); ++i) { @@ -90,18 +92,19 @@ void NormalFilter::filter(PointView& view) auto B = eigen::computeCovariance(view, ids); // perform the eigen decomposition - SelfAdjointEigenSolver solver(B); - if (solver.info() != Success) + Eigen::SelfAdjointEigenSolver solver(B); + if (solver.info() != Eigen::Success) throwError("Cannot perform eigen decomposition."); auto eval = solver.eigenvalues(); auto evec = solver.eigenvectors().col(0); - view.setField(m_nx, i, evec[0]); - view.setField(m_ny, i, evec[1]); - view.setField(m_nz, i, evec[2]); + view.setField(Dimension::Id::NormalX, i, evec[0]); + view.setField(Dimension::Id::NormalY, i, evec[1]); + view.setField(Dimension::Id::NormalZ, i, evec[2]); double sum = eval[0] + eval[1] + eval[2]; - view.setField(m_curvature, i, sum ? std::fabs(eval[0] / sum) : 0); + view.setField(Dimension::Id::Curvature, i, + sum ? std::fabs(eval[0] / sum) : 0); } } diff --git a/filters/NormalFilter.hpp b/filters/NormalFilter.hpp index 052a4a28bf..c804c02d4b 100644 --- a/filters/NormalFilter.hpp +++ b/filters/NormalFilter.hpp @@ -59,16 +59,17 @@ class PDAL_DLL NormalFilter : public Filter NormalFilter& operator=(const NormalFilter&) = delete; NormalFilter(const NormalFilter&) = delete; + void doFilter(PointView& view); + static void * create(); static int32_t destroy(void *); std::string getName() const; private: int m_knn; - Dimension::Id m_nx, m_ny, m_nz, m_curvature; - virtual void addDimensions(PointLayoutPtr layout); virtual void addArgs(ProgramArgs& args); + virtual void addDimensions(PointLayoutPtr layout); virtual void filter(PointView& view); }; diff --git a/io/PlyWriter.cpp b/io/PlyWriter.cpp index 407a8e1a05..89564276d8 100644 --- a/io/PlyWriter.cpp +++ b/io/PlyWriter.cpp @@ -34,6 +34,7 @@ #include "PlyWriter.hpp" +#include #include #include @@ -92,6 +93,8 @@ void PlyWriter::addArgs(ProgramArgs& args) { args.add("filename", "Output filename", m_filename).setPositional(); args.add("storage_mode", "PLY Storage mode", m_storageModeSpec, "default"); + args.add("dims", "Dimension names", m_dimNames); + args.add("faces", "Write faces", m_faces); } @@ -125,6 +128,28 @@ void PlyWriter::initialize() } +void PlyWriter::prepared(PointTableRef table) +{ + if (m_dimNames.size()) + { + for (auto& name : m_dimNames) + { + auto id = table.layout()->findDim(name); + if (id == Dimension::Id::Unknown) + throwError("Unknown dimension '" + name + "' in provided " + "dimension list."); + m_dims.push_back(id); + } + } + else + { + m_dims = table.layout()->dims(); + for (auto dim : m_dims) + m_dimNames.push_back(Utils::tolower(table.layout()->dimName(dim))); + } +} + + void PlyWriter::ready(PointTableRef table) { try @@ -144,7 +169,20 @@ void PlyWriter::ready(PointTableRef table) void PlyWriter::write(const PointViewPtr data) { + if (m_faces && m_pointCollector->size()) + throwError("Can't output faces for pipeline with multiple " + "point views."); m_pointCollector->append(*data); + if (m_faces) + { + m_mesh = data->mesh(); + if (!m_mesh) + throwError("Can't write mesh faces. No mesh has been generated."); + } + if (data->size() > std::numeric_limits::max()) + throwError("Can't write PLY file. Only " + + std::to_string(std::numeric_limits::max()) + + " points supported."); } @@ -155,13 +193,19 @@ void PlyWriter::done(PointTableRef table) if (!ply_add_element(m_ply, "vertex", m_pointCollector->size())) throwError("Could not add vertex element"); - auto dimensions = table.layout()->dims(); - for (auto dim : dimensions) { - std::string name = table.layout()->dimName(dim); + auto ni = m_dimNames.begin(); + for (auto dim : m_dims) { + std::string name = *ni++; e_ply_type plyType = getPlyType(table.layout()->dimType(dim)); if (!ply_add_scalar_property(m_ply, name.c_str(), plyType)) throwError("Could not add scalar property '" + name + "'"); } + if (m_faces) + { + ply_add_element(m_ply, "face", m_mesh->size()); + ply_add_list_property(m_ply, "vertex_indices", PLY_UINT8, PLY_UINT32); + } + if (!ply_add_comment(m_ply, "Generated by PDAL")) throwError("Could not add comment"); @@ -170,12 +214,26 @@ void PlyWriter::done(PointTableRef table) for (PointId index = 0; index < m_pointCollector->size(); ++index) { - for (auto dim : dimensions) + auto ni = m_dimNames.begin(); + auto di = m_dims.begin(); + for (; di != m_dims.end(); ++di, ++ni) { - double value = m_pointCollector->getFieldAs(dim, index); + double value = m_pointCollector->getFieldAs(*di, index); if (!ply_write(m_ply, value)) - throwError("Error writing dimension '" + - table.layout()->dimName(dim) + "' of point number " + + throwError("Error writing dimension '" + *ni + "' of point " + "number " + Utils::toString(index) + "."); + } + } + if (m_faces) + { + for (PointId id = 0; id < m_mesh->size(); id++) + { + const Triangle& t = (*m_mesh)[id]; + if (!(ply_write(m_ply, 3) && + ply_write(m_ply, t.m_a) && + ply_write(m_ply, t.m_b) && + ply_write(m_ply, t.m_c))) + throwError("Error writing face number " + Utils::toString(index) + "."); } } diff --git a/io/PlyWriter.hpp b/io/PlyWriter.hpp index 39dff77593..f8b8fc6a30 100644 --- a/io/PlyWriter.hpp +++ b/io/PlyWriter.hpp @@ -43,6 +43,7 @@ extern "C" PF_ExitFunc PlyWriter_InitPlugin(); namespace pdal { +class TriangularMesh; class PDAL_DLL PlyWriter : public Writer { @@ -63,6 +64,7 @@ class PDAL_DLL PlyWriter : public Writer private: virtual void addArgs(ProgramArgs& args); virtual void initialize(); + virtual void prepared(PointTableRef table); virtual void ready(PointTableRef table); virtual void write(const PointViewPtr data); virtual void done(PointTableRef table); @@ -72,6 +74,10 @@ class PDAL_DLL PlyWriter : public Writer PointViewPtr m_pointCollector; std::string m_storageModeSpec; e_ply_storage_mode m_storageMode; + bool m_faces; + StringList m_dimNames; + Dimension::IdList m_dims; + TriangularMesh *m_mesh; }; } diff --git a/pdal/Dimension.json b/pdal/Dimension.json index a44ab14e25..3d381fd34b 100644 --- a/pdal/Dimension.json +++ b/pdal/Dimension.json @@ -30,6 +30,11 @@ "description": "Z component of normal vector" }, { + "name": "Curvature", + "type": "double", + "description": "Curvature of surface at point" + }, + { "name": "Intensity", "type": "uint16", "description": "Representation of the pulse return magnitude" diff --git a/pdal/Mesh.hpp b/pdal/Mesh.hpp index afe26d228d..41e1db4af0 100644 --- a/pdal/Mesh.hpp +++ b/pdal/Mesh.hpp @@ -71,6 +71,8 @@ class PDAL_DLL TriangularMesh : public Mesh { return m_index.size(); } void add(PointId a, PointId b, PointId c) { m_index.emplace_back(a, b, c); } + const Triangle& operator[](PointId id) const + { return m_index[id]; } protected: std::deque m_index; }; diff --git a/pdal/PointView.cpp b/pdal/PointView.cpp index d3574b3ed7..19233b4081 100644 --- a/pdal/PointView.cpp +++ b/pdal/PointView.cpp @@ -34,6 +34,7 @@ #include +#include #include #include #include @@ -56,10 +57,7 @@ PointView::PointView(PointTableRef pointTable, const SpatialReference& srs) : } PointView::~PointView() -{ - for (auto it = m_meshes.begin(); it != m_meshes.end(); ++it) - delete it->second; -} +{} PointViewIter PointView::begin() { @@ -166,9 +164,10 @@ TriangularMesh *PointView::createMesh(const std::string& name) { if (Utils::contains(m_meshes, name)) return nullptr; - auto res = m_meshes.insert(std::make_pair(name, new TriangularMesh)); + auto res = m_meshes.insert(std::make_pair(name, + std::unique_ptr(new TriangularMesh))); if (res.second) - return res.first->second; + return res.first->second.get(); return nullptr; } @@ -177,11 +176,42 @@ TriangularMesh *PointView::mesh(const std::string& name) { auto it = m_meshes.find(name); if (it != m_meshes.end()) - return it->second; + return it->second.get(); + if (name.empty() && m_meshes.size()) + return m_meshes.begin()->second.get(); return nullptr; } +KD3Index& PointView::build3dIndex() +{ + //ABELL + // Should we allow a force of point view build - perhaps the index has + // changed or the point values have changed. + if (!m_index3) + { + m_index3.reset(new KD3Index(*this)); + std::cerr << "About to build index!\n"; + m_index3->build(); + } + return *m_index3.get(); +} + + +KD2Index& PointView::build2dIndex() +{ + //ABELL + // Should we allow a force of point view build - perhaps the index has + // changed or the point values have changed. + if (!m_index2) + { + m_index2.reset(new KD2Index(*this)); + m_index2->build(); + } + return *m_index2.get(); +} + + void PointView::dump(std::ostream& ostr) const { using std::endl; diff --git a/pdal/PointView.hpp b/pdal/PointView.hpp index 1d87222f30..dce1ceded6 100644 --- a/pdal/PointView.hpp +++ b/pdal/PointView.hpp @@ -63,6 +63,8 @@ namespace plang struct PointViewLess; class PointView; class PointViewIter; +class KD2Index; +class KD3Index; typedef std::shared_ptr PointViewPtr; typedef std::set PointViewSet; @@ -73,6 +75,8 @@ class PDAL_DLL PointView : public PointContainer friend class PointIdxRef; friend struct PointViewLess; public: + PointView(const PointView&) = delete; + PointView& operator=(const PointView&) = delete; PointView(PointTableRef pointTable); PointView(PointTableRef pointTable, const SpatialReference& srs); @@ -243,7 +247,6 @@ class PDAL_DLL PointView : public PointContainer } } - /// Provides access to the memory storing the point data. Though this /// function is public, other access methods are safer and preferred. char *getPoint(PointId id) @@ -288,7 +291,10 @@ class PDAL_DLL PointView : public PointContainer \param name Name of the mesh. \return New mesh. Null is returned if the mesh already exists. */ - TriangularMesh *mesh(const std::string& name); + TriangularMesh *mesh(const std::string& name = ""); + + KD3Index& build3dIndex(); + KD2Index& build2dIndex(); protected: PointTableRef m_pointTable; @@ -299,7 +305,9 @@ class PDAL_DLL PointView : public PointContainer int m_id; std::queue m_temps; SpatialReference m_spatialReference; - std::map m_meshes; + std::map> m_meshes; + std::unique_ptr m_index3; + std::unique_ptr m_index2; private: static int m_lastId; diff --git a/test/unit/PointViewTest.cpp b/test/unit/PointViewTest.cpp index 16b9359ebf..f9b9eb8504 100644 --- a/test/unit/PointViewTest.cpp +++ b/test/unit/PointViewTest.cpp @@ -161,58 +161,6 @@ TEST(PointViewTest, getFloat) } -TEST(PointViewTest, copy) -{ - PointTable table; - PointViewPtr view = makeTestView(table); - - PointView d2(*view); - - // read the view back out - { - EXPECT_EQ( - d2.getFieldAs(Dimension::Id::Classification, 0), - view->getFieldAs(Dimension::Id::Classification, 0)); - EXPECT_EQ(d2.getFieldAs(Dimension::Id::X, 0), - view->getFieldAs(Dimension::Id::X, 0)); - EXPECT_FLOAT_EQ(d2.getFieldAs(Dimension::Id::Y, 0), - view->getFieldAs(Dimension::Id::Y, 0)); - } - - for (int i = 1; i < 17; i++) - { - uint8_t x = d2.getFieldAs(Dimension::Id::Classification, i); - int32_t y = d2.getFieldAs(Dimension::Id::X, i); - double z = d2.getFieldAs(Dimension::Id::Y, i); - - EXPECT_EQ(x, i + 1u); - EXPECT_EQ(y, i * 10); - EXPECT_TRUE(Utils::compare_approx(z, i * 100.0, - (std::numeric_limits::min)())); - } - EXPECT_EQ(view->size(), d2.size()); - -} - -TEST(PointViewTest, copyCtor) -{ - PointTable table; - PointViewPtr view = makeTestView(table); - - PointView d2(*view); - verifyTestView(d2); -} - -TEST(PointViewTest, assignment) -{ - PointTable table; - PointViewPtr view = makeTestView(table); - - PointView d2 = *view; - verifyTestView(d2); -} - - TEST(PointViewTest, bigfile) { PointTable table; @@ -280,68 +228,6 @@ TEST(PointViewTest, bigfile) } -//ABELL - Move to KdIndex -/** -TEST(PointViewTest, kdindex) -{ - LasReader reader(Support::viewpath("1.2-with-color.las")); - reader.prepare(); - - const Schema& schema = reader.getSchema(); - uint32_t capacity(1000); - PointView view(schema, capacity); - - StageSequentialIterator* iter = reader.createSequentialIterator(view); - - { - uint32_t numRead = iter->read(view); - EXPECT_EQ(numRead, capacity); - } - - EXPECT_EQ(view.getCapacity(), capacity); - EXPECT_EQ(view.getSchema(), schema); - - - IndexedPointView iview(view); - EXPECT_EQ(iview.getCapacity(), capacity); - EXPECT_EQ(iview.getSchema(), schema); - - iview.build(); - - unsigned k = 8; - - // If the query distance is 0, just return the k nearest neighbors - std::vector ids = idata.neighbors(636199, 849238, 428.05, k); - EXPECT_EQ(ids.size(), k); - EXPECT_EQ(ids[0], 8u); - EXPECT_EQ(ids[1], 7u); - EXPECT_EQ(ids[2], 9u); - EXPECT_EQ(ids[3], 42u); - EXPECT_EQ(ids[4], 40u); - - std::vector dist_ids = idata.neighbors(636199, 849238, 428.05, 3); - - EXPECT_EQ(dist_ids.size(), 3u); - EXPECT_EQ(dist_ids[0], 8u); - - std::vector nids = idata.neighbors(636199, 849238, 428.05, k); - - EXPECT_EQ(nids.size(), k); - EXPECT_EQ(nids[0], 8u); - EXPECT_EQ(nids[1], 7u); - EXPECT_EQ(nids[2], 9u); - EXPECT_EQ(nids[3], 42u); - EXPECT_EQ(nids[4], 40u); - - std::vector rids = iview.radius(637012.24, 849028.31, - 431.66, 100000); - EXPECT_EQ(rids.size(), 11u); - - delete iter; -} -**/ - - static void check_bounds(const BOX3D& box, double minx, double maxx, double miny, double maxy, From d03984916e8ebdba62ae5d3b1688a93cb7228d83 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Thu, 20 Apr 2017 12:30:54 -0500 Subject: [PATCH 20/57] Add mesh writer test. --- test/data/ply/mesh.ply | 16 +++++++++++ test/unit/io/PlyWriterTest.cpp | 52 ++++++++++++++++++++++++++++++++++ 2 files changed, 68 insertions(+) create mode 100644 test/data/ply/mesh.ply diff --git a/test/data/ply/mesh.ply b/test/data/ply/mesh.ply new file mode 100644 index 0000000000..42c921af5e --- /dev/null +++ b/test/data/ply/mesh.ply @@ -0,0 +1,16 @@ +ply +format ascii 1.0 +comment Generated by PDAL +element vertex 4 +property float64 x +property float64 y +property float64 z +element face 2 +property list uint8 uint32 vertex_indices +end_header +1 1 0 +2 1 0 +1 2 0 +2 2 2 +3 0 1 2 +3 1 2 3 diff --git a/test/unit/io/PlyWriterTest.cpp b/test/unit/io/PlyWriterTest.cpp index 8ed588ff9f..dd4acf2cad 100644 --- a/test/unit/io/PlyWriterTest.cpp +++ b/test/unit/io/PlyWriterTest.cpp @@ -35,6 +35,8 @@ #include #include +#include +#include #include #include #include "Support.hpp" @@ -73,5 +75,55 @@ TEST(PlyWriter, Write) writer.execute(table); } +TEST(PlyWriter, mesh) +{ + std::string outfile(Support::temppath("out.ply")); + std::string testfile(Support::datapath("ply/mesh.ply")); + + FileUtils::deleteFile(outfile); + + PointTable t; + + t.layout()->registerDim(Dimension::Id::X); + t.layout()->registerDim(Dimension::Id::Y); + t.layout()->registerDim(Dimension::Id::Z); + + PointViewPtr v(new PointView(t)); + v->setField(Dimension::Id::X, 0, 1); + v->setField(Dimension::Id::Y, 0, 1); + v->setField(Dimension::Id::Z, 0, 0); + + v->setField(Dimension::Id::X, 1, 2); + v->setField(Dimension::Id::Y, 1, 1); + v->setField(Dimension::Id::Z, 1, 0); + + v->setField(Dimension::Id::X, 2, 1); + v->setField(Dimension::Id::Y, 2, 2); + v->setField(Dimension::Id::Z, 2, 0); + + v->setField(Dimension::Id::X, 3, 2); + v->setField(Dimension::Id::Y, 3, 2); + v->setField(Dimension::Id::Z, 3, 2); + + TriangularMesh *mesh = v->createMesh("foo"); + mesh->add(0, 1, 2); + mesh->add(1, 2, 3); + + BufferReader r; + r.addView(v); + + PlyWriter w; + Options wo; + wo.add("filename", outfile); + wo.add("storage_mode", "ascii"); + wo.add("faces", true); + w.setInput(r); + w.setOptions(wo); + + w.prepare(t); + w.execute(t); + + EXPECT_TRUE(Support::compare_text_files(outfile, testfile)); +} } From d939209b9e8a4a481a53bef64ff4ecfb3bf89b38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Lemoine?= Date: Fri, 21 Apr 2017 14:37:31 +0200 Subject: [PATCH 21/57] Fix formatting issues in testing.rst (#1574) --- doc/development/testing.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/development/testing.rst b/doc/development/testing.rst index 7997a1cffb..ac8f3cdd42 100644 --- a/doc/development/testing.rst +++ b/doc/development/testing.rst @@ -315,14 +315,14 @@ and to run only select tests (``--gtest_filter``). .. note:: If the PostgreSQL PointCloud plugin was enabled on the CMake command line - (with ``-DBUILD_PLUGIN_PGPOINTCLOUD=ON``) then `ctest` will attempt to run - the ``pgpointcloud`` tests. And you will get PostgreSQL connection errors + (with ``-DBUILD_PLUGIN_PGPOINTCLOUD=ON``) then ``ctest`` will attempt to run + the ``pgpointcloud`` tests. And you will get PostgreSQL connection errors if the `libpq environment variables`_ are not correctly set in your shell. This is for example how you can run the ``pgpointcloud`` tests:: $ PGUSER=pdal PGPASSWORD=pdal PGHOST=localhost ctest -R pgpointcloudtest -.. _`libpg environment variables`: https://www.postgresql.org/docs/current/static/libpq-envars.html +.. _`libpq environment variables`: https://www.postgresql.org/docs/current/static/libpq-envars.html Test Data ========= From 92184d076ed1a0d25c5402f0a0e0bf2a5e451333 Mon Sep 17 00:00:00 2001 From: Howard Butler Date: Fri, 21 Apr 2017 16:28:56 -0500 Subject: [PATCH 22/57] we need to count clusters from 1 instead of 0, because all other points will have 0 for their ClusterID (default value) (#1576) --- filters/ClusterFilter.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/filters/ClusterFilter.cpp b/filters/ClusterFilter.cpp index ffe1aa8281..0ac35a367d 100644 --- a/filters/ClusterFilter.cpp +++ b/filters/ClusterFilter.cpp @@ -73,7 +73,7 @@ void ClusterFilter::filter(PointView& view) auto clusters = Segmentation::extractClusters(view, m_minPoints, m_maxPoints, m_tolerance); - uint64_t id = 0; + uint64_t id = 1; for (auto const& c : clusters) { for (auto const& i : c) From e6c0f8fe4ec72f0bf8a13616c16372635aaca0c5 Mon Sep 17 00:00:00 2001 From: Howard Butler Date: Mon, 24 Apr 2017 15:13:28 -0500 Subject: [PATCH 23/57] don't purge python-numpy in diet section --- scripts/docker/Dockerfile | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/docker/Dockerfile b/scripts/docker/Dockerfile index 8949aba323..f241cd9f51 100644 --- a/scripts/docker/Dockerfile +++ b/scripts/docker/Dockerfile @@ -78,7 +78,6 @@ RUN apt-get purge -y \ libtiff5-dev \ openssh-client \ python-dev \ - python-numpy \ python-software-properties \ software-properties-common \ wget \ From 5889e855965b0eb312363e91cf4159d969295284 Mon Sep 17 00:00:00 2001 From: Howard Butler Date: Tue, 25 Apr 2017 09:54:11 -0500 Subject: [PATCH 24/57] more python numpy docker fixes --- scripts/docker/Dockerfile | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/scripts/docker/Dockerfile b/scripts/docker/Dockerfile index f241cd9f51..6767b04ad8 100644 --- a/scripts/docker/Dockerfile +++ b/scripts/docker/Dockerfile @@ -6,6 +6,8 @@ ENV CXX g++ RUN apt-get update && apt-get install -y --fix-missing --no-install-recommends \ libhpdf-dev \ + python-all-dev \ + python-numpy \ libsqlite3-mod-spatialite \ && rm -rf /var/lib/apt/lists/* @@ -77,7 +79,6 @@ RUN apt-get purge -y \ liblapack-dev \ libtiff5-dev \ openssh-client \ - python-dev \ python-software-properties \ software-properties-common \ wget \ @@ -114,8 +115,7 @@ RUN apt-get purge -y \ libcurl4-openssl-dev \ libspatialite-dev \ libdap-dev\ - cython \ - python-pip + cython RUN apt-get autoremove -y @@ -172,5 +172,7 @@ RUN apt-get update && apt-get install -y \ libpcl-segmentation1.7 \ libpcl-surface1.7 \ libpcl-tracking1.7 \ - libpcl-visualization1.7 + libpcl-visualization1.7 \ + python-all-dev \ + python-pip From 8de8c4c721604194ecc2198f23c632f9c9dabd2f Mon Sep 17 00:00:00 2001 From: Connor Manning Date: Tue, 25 Apr 2017 16:10:28 -0400 Subject: [PATCH 25/57] Map .csv files to the text reader. Tweak text reader for more flexibility in custom layouts. --- io/TextReader.cpp | 3 ++- pdal/StageFactory.cpp | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/io/TextReader.cpp b/io/TextReader.cpp index 40122a8f57..eaa9381247 100644 --- a/io/TextReader.cpp +++ b/io/TextReader.cpp @@ -94,12 +94,13 @@ void TextReader::addArgs(ProgramArgs& args) void TextReader::addDimensions(PointLayoutPtr layout) { + m_dims.clear(); for (auto name : m_dimNames) { Utils::trim(name); Dimension::Id id = layout->registerOrAssignDim(name, Dimension::Type::Double); - if (Utils::contains(m_dims, id)) + if (Utils::contains(m_dims, id) && id != pdal::Dimension::Id::Unknown) throwError("Duplicate dimension '" + name + "' detected in input file '" + m_filename + "'."); m_dims.push_back(id); diff --git a/pdal/StageFactory.cpp b/pdal/StageFactory.cpp index 676df532ee..42116fd365 100644 --- a/pdal/StageFactory.cpp +++ b/pdal/StageFactory.cpp @@ -126,7 +126,7 @@ StringList StageFactory::extensions(const std::string& driver) { "readers.sqlite", { "sqlite" } }, { "readers.mrsid", { "sid" } }, { "readers.tindex", { "tindex" } }, - { "readers.text", { "txt" } }, + { "readers.text", { "csv", "txt" } }, { "readers.icebridge", { "h5" } }, { "writers.bpf", { "bpf" } }, @@ -152,6 +152,7 @@ std::string StageFactory::inferReaderDriver(const std::string& filename) { "bin", "readers.terrasolid" }, { "bpf", "readers.bpf" }, { "csd", "readers.optech" }, + { "csv", "readers.text" }, { "greyhound", "readers.greyhound" }, { "icebridge", "readers.icebridge" }, { "las", "readers.las" }, From 180facd20079877f7f97b06f2711d7fa26bb45cb Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Tue, 2 May 2017 08:10:55 -0500 Subject: [PATCH 26/57] PLY reader without rply. --- io/PlyReader.cpp | 441 ++++++++++++++++++++++++++++-------------- io/PlyReader.hpp | 99 ++++++++-- pdal/PDALUtils.hpp | 5 +- pdal/PointRef.hpp | 2 + pdal/util/IStream.hpp | 163 +++++++++++++++- pdal/util/Utils.hpp | 4 +- 6 files changed, 552 insertions(+), 162 deletions(-) diff --git a/io/PlyReader.cpp b/io/PlyReader.cpp index 7f5e78a936..123c04fed0 100644 --- a/io/PlyReader.cpp +++ b/io/PlyReader.cpp @@ -36,75 +36,231 @@ #include +#include #include #include +#include namespace pdal { -namespace + +std::string PlyReader::readLine() { + m_line.clear(); + if (m_lines.size()) + { + m_line = m_lines.top(); + m_lines.pop(); + } + else + { + do + { + std::getline(*m_stream, m_line); + } while (m_line.empty() && m_stream->good()); + } + Utils::trimTrailing(m_line); + m_linePos = Utils::extract(m_line, 0, + [](char c){ return !std::isspace(c); }); + return std::string(m_line, 0, m_linePos); +} + + +void PlyReader::pushLine() +{ + m_lines.push(m_line); +} -struct CallbackContext + +std::string PlyReader::nextWord() { - PointViewPtr view; - PlyReader::DimensionMap dimensionMap; -}; + std::string s; + std::string::size_type cnt = + Utils::extract(m_line, m_linePos, (int(*)(int))std::isspace); + m_linePos += cnt; + if (m_linePos == m_line.size()) + return s; + + cnt = Utils::extract(m_line, m_linePos, + [](char c){ return !std::isspace(c); }); + s = std::string(m_line, m_linePos, cnt); + m_linePos += cnt; + return s; +} -void plyErrorCallback(p_ply ply, const char * message) +void PlyReader::extractMagic() { - throw PlyReader::error(message); + std::string first = readLine(); + if (first != "ply") + throwError("File isn't a PLY file. 'ply' not found."); + if (m_linePos != m_line.size()) + throwError("Text found following 'ply' keyword."); } -p_ply openPly(std::string filename) +void PlyReader::extractEnd() { - p_ply ply = ply_open(filename.c_str(), &plyErrorCallback, 0, nullptr); - if (!ply) - throw PlyReader::error("Unable to open file " + filename + - " for reading."); - - if (!ply_read_header(ply)) - throw PlyReader::error("Unable to read header of " + filename + "."); - return ply; + std::string first = readLine(); + if (first != "end_header") + throwError("'end_header' expected but found line beginning with '" + + first + "' instead."); + if (m_linePos != m_line.size()) + throwError("Text found following 'end_header' keyword."); } -int readPlyCallback(p_ply_argument argument) +void PlyReader::extractFormat() +{ + std::string word = readLine(); + if (word != "format") + throwError("Expected format line not found in PLY file."); + + word = nextWord(); + if (word == "ascii") + m_format = Format::ASCII; + else if (word == "binary_big_endian") + m_format = Format::BINARY_BE; + else if (word == "binary_little_endian") + m_format = Format::BINARY_LE; + else + throwError("Unrecognized PLY format: '" + word + "'."); + + word = nextWord(); + if (word != "1.0") + throwError("Unsupported PLY version: '" + word + "'."); +} + +Dimension::Type PlyReader::getType(const std::string& name) { - p_ply_element element; - long index; - void * contextAsVoid; - long numToRead; - p_ply_property property; - const char * propertyName; + static std::map types = + { + { "int8", Dimension::Type::Signed8 }, + { "uint8", Dimension::Type::Unsigned8 }, + { "int16", Dimension::Type::Signed16 }, + { "uint16", Dimension::Type::Unsigned16 }, + { "int32", Dimension::Type::Signed32 }, + { "uint32", Dimension::Type::Unsigned32 }, + { "float32", Dimension::Type::Float }, + { "float64", Dimension::Type::Double }, + + { "char", Dimension::Type::Signed8 }, + { "uchar", Dimension::Type::Unsigned8 }, + { "short", Dimension::Type::Signed16 }, + { "ushort", Dimension::Type::Unsigned16 }, + { "int", Dimension::Type::Signed32 }, + { "uint", Dimension::Type::Unsigned32 }, + { "float", Dimension::Type::Float }, + { "double", Dimension::Type::Double } + }; + + try + { + return types.at(name); + } + catch (std::out_of_range) + {} + return Dimension::Type::None; +} - if (!ply_get_argument_element(argument, &element, &index)) - throw PlyReader::error("Error getting argument element."); - if (!ply_get_argument_user_data(argument, &contextAsVoid, &numToRead)) - throw PlyReader::error("Error getting argument user data."); +void PlyReader::extractProperty(Element& element) +{ + std::string word = nextWord(); + Dimension::Type type = getType(word); - // We've read enough, abort the callback cycle - if (numToRead <= index) - return 0; + if (type != Dimension::Type::None) + { + std::string name = nextWord(); + if (name.empty()) + throwError("No name for property of element '" + + element.m_name + "'."); + element.m_properties.push_back( + std::unique_ptr(new SimpleProperty(name, type))); + } + else if (word == "list") + { + if (element.m_name == "vertex") + throwError("List properties are not supported for the 'vertex' " + "element."); + + word = nextWord(); + Dimension::Type countType = getType(word); + if (countType == Dimension::Type::None) + throwError("No valid count type for list property of element '" + + element.m_name + "'."); + word = nextWord(); + Dimension::Type listType = getType(word); + if (listType == Dimension::Type::None) + throwError("No valid list type for list property of element '" + + element.m_name + "'."); + std::string name = nextWord(); + if (name.empty()) + throwError("No name for property of element '" + + element.m_name + "'."); + element.m_properties.push_back( + std::unique_ptr(new ListProperty(name, countType, + listType))); + } + else + throwError("Invalid property type '" + word + "'."); +} - if (!ply_get_argument_property(argument, &property, nullptr, nullptr)) - throw PlyReader::error("Error getting argument property."); - if (!ply_get_property_info(property, &propertyName, nullptr, - nullptr, nullptr)) - throw PlyReader::error("Error getting property info."); +void PlyReader::extractProperties(Element& element) +{ + while (true) + { + std::string word = readLine(); + if (word == "comment" || word == "obj_info") + continue; + else if (word == "property") + extractProperty(element); + else + { + pushLine(); + break; + } + } +} - CallbackContext * context = static_cast(contextAsVoid); - double value = ply_get_argument_value(argument); - Dimension::Id dimension = context->dimensionMap.at(propertyName); - context->view->setField(dimension, index, value); - return 1; +bool PlyReader::extractElement() +{ + std::string word = readLine(); + if (word == "comment" || word == "obj_info") + return true; + else if (word == "end_header") + { + pushLine(); + return false; + } + else if (word == "element") + { + std::string name = nextWord(); + if (name.empty()) + throwError("Missing element name."); + long count = std::stol(nextWord()); + if (count < 0) + throwError("Invalid count for element '" + name + "'."); + m_elements.emplace_back(name, count); + extractProperties(m_elements.back()); + return true; + } + throwError("Invalid keyword '" + word + "' when expecting an element."); + return false; // quiet compiler } -} // unnamed namespace + +void PlyReader::extractHeader() +{ + extractMagic(); + extractFormat(); + while (extractElement()) + ; + extractEnd(); + m_dataPos = m_stream->tellg(); +} static PluginInfo const s_info = PluginInfo( @@ -123,147 +279,148 @@ std::string PlyReader::getName() const PlyReader::PlyReader() - : m_ply(nullptr) - , m_vertexDimensions() {} void PlyReader::initialize() { - try + m_stream = Utils::openFile(m_filename, true); + extractHeader(); + Utils::closeFile(m_stream); + m_stream = nullptr; +} + + +void PlyReader::addDimensions(PointLayoutPtr layout) +{ + SimpleProperty *prop; + + for (auto& elt : m_elements) { - p_ply ply = openPly(m_filename); - p_ply_element vertex_element = nullptr; - bool found_vertex_element = false; - const char* element_name; - long element_count; - while ((vertex_element = ply_get_next_element(ply, vertex_element))) + if (elt.m_name == "vertex") { - if (!ply_get_element_info(vertex_element, &element_name, - &element_count)) - throwError("Error reading element info in " + m_filename + "."); - if (strcmp(element_name, "vertex") == 0) + for (auto& prop : elt.m_properties) { - found_vertex_element = true; - break; + auto vprop = static_cast(prop.get()); + layout->registerOrAssignDim(vprop->m_name, vprop->m_type); + vprop->setDim( + layout->registerOrAssignDim(vprop->m_name, vprop->m_type)); } + return; } - if (!found_vertex_element) - throwError("File " + m_filename + " does not contain a vertex " - "element."); + } + throwError("No 'vertex' element in header."); +} - static std::map types = - { - { PLY_INT8, Dimension::Type::Signed8 }, - { PLY_UINT8, Dimension::Type::Unsigned8 }, - { PLY_INT16, Dimension::Type::Signed16 }, - { PLY_UINT16, Dimension::Type::Unsigned16 }, - { PLY_INT32, Dimension::Type::Signed32 }, - { PLY_UINT32, Dimension::Type::Unsigned32 }, - { PLY_FLOAT32, Dimension::Type::Float }, - { PLY_FLOAT64, Dimension::Type::Double }, - - { PLY_CHAR, Dimension::Type::Signed8 }, - { PLY_UCHAR, Dimension::Type::Unsigned8 }, - { PLY_SHORT, Dimension::Type::Signed16 }, - { PLY_USHORT, Dimension::Type::Unsigned16 }, - { PLY_INT, Dimension::Type::Signed32 }, - { PLY_UINT, Dimension::Type::Unsigned32 }, - { PLY_FLOAT, Dimension::Type::Float }, - { PLY_DOUBLE, Dimension::Type::Double } - }; - - p_ply_property property = nullptr; - while ((property = ply_get_next_property(vertex_element, property))) - { - const char* name; - e_ply_type type; - e_ply_type length_type; - e_ply_type value_type; - if (!ply_get_property_info(property, &name, &type, - &length_type, &value_type)) - throwError("Error reading property info in " + - m_filename + "."); - // For now, we'll just use PDAL's built in dimension matching. - // We could be smarter about this, e.g. by using the length - // and value type attributes. - m_vertexTypes[name] = types[type]; - } - ply_close(ply); + +bool PlyReader::readProperty(Property *prop, PointRef& point) +{ + if (!m_stream->good()) + return false; + prop->read(m_stream, m_format, point); + return true; +} + + +void PlyReader::SimpleProperty::read(std::istream *stream, + PlyReader::Format format, PointRef& point) +{ + if (format == Format::ASCII) + { + double d; + *stream >> d; + point.setField(m_dim, d); } - catch (const error& err) + else if (format == Format::BINARY_LE) { - throwError(err.what()); + ILeStream in(stream); + Everything e = Utils::extractDim(in, m_type); + point.setField(m_dim, m_type, &e); + } + else if (format == Format::BINARY_BE) + { + IBeStream in(stream); + Everything e = Utils::extractDim(in, m_type); + point.setField(m_dim, m_type, &e); } } -void PlyReader::addDimensions(PointLayoutPtr layout) +// Right now we don't support list properties for point data. We just +// read the data and throw it away. +void PlyReader::ListProperty::read(std::istream *stream, + PlyReader::Format format, PointRef& point) { - for (auto it : m_vertexTypes) + if (format == Format::ASCII) { - const std::string& name = it.first; - const Dimension::Type& type = it.second; + size_t cnt; + *stream >> cnt; - m_vertexDimensions[name] = layout->registerOrAssignDim(name, type); + double d; + while (cnt--) + *stream >> d; + } + else if (format == Format::BINARY_LE) + { + ILeStream istream(stream); + Everything e = Utils::extractDim(istream, m_countType); + size_t cnt = (size_t)Utils::toDouble(e, m_countType); + cnt *= Dimension::size(m_listType); + istream.seek(cnt, std::ios_base::cur); } + else if (format == Format::BINARY_BE) + { + IBeStream istream(stream); + Everything e = Utils::extractDim(istream, m_countType); + size_t cnt = (size_t)Utils::toDouble(e, m_countType); + cnt *= Dimension::size(m_listType); + istream.seek(cnt, std::ios_base::cur); + } +} + + +void PlyReader::readElement(Element& elt, PointRef& point) +{ + for (auto& prop : elt.m_properties) + if (!readProperty(prop.get(), point)) + throwError("Error reading data for point " + + std::to_string(point.pointId()) + "."); } void PlyReader::ready(PointTableRef table) { - try - { - m_ply = openPly(m_filename); - } - catch (const error& err) - { - throwError(err.what()); - } + m_stream = Utils::openFile(m_filename, true); + if (m_stream) + m_stream->seekg(m_dataPos); } point_count_t PlyReader::read(PointViewPtr view, point_count_t num) { - CallbackContext context; - context.view = view; - context.dimensionMap = m_vertexDimensions; - - // It's possible that point_count_t holds a value that's larger than the - // long that is the maximum rply (don't know about ply) point count. - long cnt; - cnt = Utils::inRange(num) ? num : (std::numeric_limits::max)(); - try + point_count_t cnt; + + PointRef point(view->point(0)); + for (Element& elt : m_elements) { - for (auto it : m_vertexDimensions) - { - ply_set_read_cb(m_ply, "vertex", it.first.c_str(), readPlyCallback, - &context, cnt); - } - if (!ply_read(m_ply)) + cnt = 0; + for (PointId idx = 0; idx < elt.m_count && idx < num; ++idx) { - throwError("Error reading " + m_filename + "."); + point.setPointId(idx); + readElement(elt, point); + cnt++; } + // We currently only load the vertex data. + if (elt.m_name == "vertex") + break; } - catch(const error& err) - { - throwError(err.what()); - } - return view->size(); + return cnt; } void PlyReader::done(PointTableRef table) { - try - { - if (!ply_close(m_ply)) - throwError("Error closing " + m_filename + "."); - } - catch (const error& err) - { - throwError(err.what()); - } + Utils::closeFile(m_stream); } } // namespace pdal diff --git a/io/PlyReader.hpp b/io/PlyReader.hpp index 728646e49a..3b4d8657e1 100644 --- a/io/PlyReader.hpp +++ b/io/PlyReader.hpp @@ -35,8 +35,7 @@ #pragma once #include - -#include +#include #include #include @@ -52,12 +51,6 @@ namespace pdal class PDAL_DLL PlyReader : public Reader { public: - struct error : public std::runtime_error - { - error(const std::string& e) : std::runtime_error(e) - {} - }; - static void *create(); static int32_t destroy(void *); std::string getName() const; @@ -66,19 +59,95 @@ class PDAL_DLL PlyReader : public Reader PlyReader(); - static Dimension::IdList getDefaultDimensions(); - private: + enum class Format + { + ASCII, + BINARY_LE, + BINARY_BE + }; + + struct Property + { + Property(const std::string& name) : m_name(name) + {} + + std::string m_name; + + virtual void setDim(Dimension::Id id) + {} + virtual void read(std::istream *stream, PlyReader::Format format, + PointRef& point) = 0; + }; + + struct SimpleProperty : public Property + { + SimpleProperty(const std::string& name, Dimension::Type type) : + Property(name), m_type(type), m_dim(Dimension::Id::Unknown) + {} + + Dimension::Type m_type; + Dimension::Id m_dim; + + virtual void read(std::istream *stream, PlyReader::Format format, + PointRef& point) override; + virtual void setDim(Dimension::Id id) override + { m_dim = id; } + }; + + struct ListProperty : public Property + { + ListProperty(const std::string& name, Dimension::Type countType, + Dimension::Type listType) : Property(name), m_countType(countType), + m_listType(listType) + {} + + Dimension::Type m_countType; + Dimension::Type m_listType; + + virtual void read(std::istream *stream, PlyReader::Format format, + PointRef& point) override; + }; + + struct Element + { + Element(const std::string name, size_t count) : + m_name(name), m_count(count) + {} + + std::string m_name; + size_t m_count; + std::vector> m_properties; + }; + + Format m_format; + std::string m_line; + std::string::size_type m_linePos; + std::stack m_lines; + std::istream *m_stream; + std::istream::streampos m_dataPos; + std::vector m_elements; + virtual void initialize(); virtual void addDimensions(PointLayoutPtr layout); virtual void ready(PointTableRef table); virtual point_count_t read(PointViewPtr view, point_count_t num); virtual void done(PointTableRef table); - p_ply m_ply; - - DimensionMap m_vertexDimensions; - std::map m_vertexTypes; + std::string readLine(); + void pushLine(); + std::string nextWord(); + void extractMagic(); + void extractEnd(); + void extractFormat(); + Dimension::Type getType(const std::string& name); + void extractProperty(Element& element); + void extractProperties(Element& element); + bool extractElement(); + void extractHeader(); + void readElement(Element& elt, PointRef& point); + bool readProperty(Property *prop, PointRef& point); }; -} + +} // namespace pdal diff --git a/pdal/PDALUtils.hpp b/pdal/PDALUtils.hpp index 0d866335d1..90af65d1bd 100644 --- a/pdal/PDALUtils.hpp +++ b/pdal/PDALUtils.hpp @@ -110,7 +110,8 @@ inline double toDouble(const Everything& e, Dimension::Type type) return d; } -inline Everything extractDim(Extractor& ext, Dimension::Type type) +template +inline Everything extractDim(IN& ext, Dimension::Type type) { using Type = Dimension::Type; @@ -153,6 +154,7 @@ inline Everything extractDim(Extractor& ext, Dimension::Type type) return e; } + inline void insertDim(Inserter& ins, Dimension::Type type, const Everything& e) { @@ -196,7 +198,6 @@ inline void insertDim(Inserter& ins, Dimension::Type type, } - inline MetadataNode toMetadata(const BOX2D& bounds) { MetadataNode output("bbox"); diff --git a/pdal/PointRef.hpp b/pdal/PointRef.hpp index 0a201fc35e..e688c15b55 100644 --- a/pdal/PointRef.hpp +++ b/pdal/PointRef.hpp @@ -158,6 +158,8 @@ class PDAL_DLL PointRef void setPointId(PointId idx) { m_idx = idx; } + PointId pointId() const + { return m_idx; } inline void getField(char *val, Dimension::Id d, Dimension::Type type) const; inline void setField(Dimension::Id dim, diff --git a/pdal/util/IStream.hpp b/pdal/util/IStream.hpp index dd819eb7c5..1d921da302 100644 --- a/pdal/util/IStream.hpp +++ b/pdal/util/IStream.hpp @@ -420,6 +420,167 @@ class ILeStream : public IStream }; +/** + Stream wrapper for input of binary data that converts from big-endian + to host ordering. +*/ +class IBeStream : public IStream +{ +public: + /** + Default constructor. + */ + PDAL_DLL IBeStream() + {} + + /** + Constructor that opens the file and maps it to a stream. + + \param filename Filename. + */ + PDAL_DLL IBeStream(const std::string& filename) : IStream(filename) + {} + + /** + Constructor that maps to a provided stream. + + \param stream Stream to extract from. + */ + PDAL_DLL IBeStream(std::istream *stream) : IStream(stream) + {} + + /** + Extract an unsigned byte from the stream. + + \param v unsigned byte to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (uint8_t& v) + { + v = (uint8_t)m_stream->get(); + return *this; + } + + /** + Extract an unsigned byte from the stream. + + \param v unsigned byte to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (int8_t& v) + { + v = (int8_t)m_stream->get(); + return *this; + } + + /** + Extract an unsigned short from the stream. + + \param v unsigned short to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (uint16_t& v) + { + m_stream->read((char *)&v, sizeof(v)); + v = be16toh(v); + return *this; + } + + /** + Extract an short from the stream. + + \param v short to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (int16_t& v) + { + m_stream->read((char *)&v, sizeof(v)); + v = (int16_t)be16toh((uint16_t)v); + return *this; + } + + /** + Extract an unsigned int from the stream. + + \param v unsigned int to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (uint32_t& v) + { + m_stream->read((char *)&v, sizeof(v)); + v = be32toh(v); + return *this; + } + + /** + Extract an int from the stream. + + \param v int to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (int32_t& v) + { + m_stream->read((char *)&v, sizeof(v)); + v = (int32_t)be32toh((uint32_t)v); + return *this; + } + + /** + Extract an unsigned long int from the stream. + + \param v unsigned long int to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (uint64_t& v) + { + m_stream->read((char *)&v, sizeof(v)); + v = be64toh(v); + return *this; + } + + /** + Extract a long int from the stream. + + \param v long int to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (int64_t& v) + { + m_stream->read((char *)&v, sizeof(v)); + v = (int64_t)be64toh((uint64_t)v); + return *this; + } + + /** + Extract a float from the stream. + + \param v float to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (float& v) + { + m_stream->read((char *)&v, sizeof(v)); + uint32_t tmp = be32toh(*(uint32_t *)(&v)); + std::memcpy(&v, &tmp, sizeof(tmp)); + return *this; + } + + /** + Extract a double from the stream. + + \param v double to populate + \return This stream. + */ + PDAL_DLL IBeStream& operator >> (double& v) + { + m_stream->read((char *)&v, sizeof(v)); + uint64_t tmp = be64toh(*(uint64_t *)(&v)); + std::memcpy(&v, &tmp, sizeof(tmp)); + return *this; + } +}; + + /** Stream wrapper for input of binary data that converts from either little-endian or big-endian to host ordering, depending on object @@ -437,7 +598,7 @@ class ISwitchableStream : public IStream : IStream(filename) , m_isLittleEndian(DefaultIsLittleEndian) {} - + PDAL_DLL ISwitchableStream(std::istream* stream) : IStream(stream) , m_isLittleEndian(DefaultIsLittleEndian) diff --git a/pdal/util/Utils.hpp b/pdal/util/Utils.hpp index 9f5eda4d7e..0ef277392c 100644 --- a/pdal/util/Utils.hpp +++ b/pdal/util/Utils.hpp @@ -443,7 +443,7 @@ namespace Utils \param s String in which to start counting characters. \param p Position in input string at which to start counting. - \param pred Unary predicte that tests a character. + \param pred Unary predicate that tests a character. \return Then number of characters matching the predicate. */ template @@ -451,7 +451,7 @@ namespace Utils extract(const std::string& s, std::string::size_type p, PREDICATE pred) { std::string::size_type count = 0; - while (pred(s[p++])) + while (p < s.size() && pred(s[p++])) count++; return count; } From ab94bd3670652f190eba53edb04999a237d45c7c Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Tue, 2 May 2017 12:09:58 -0500 Subject: [PATCH 27/57] Change text on PLY error. --- io/PlyReader.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/io/PlyReader.cpp b/io/PlyReader.cpp index 123c04fed0..9df86df1e8 100644 --- a/io/PlyReader.cpp +++ b/io/PlyReader.cpp @@ -383,7 +383,7 @@ void PlyReader::readElement(Element& elt, PointRef& point) { for (auto& prop : elt.m_properties) if (!readProperty(prop.get(), point)) - throwError("Error reading data for point " + + throwError("Error reading data for point/element " + std::to_string(point.pointId()) + "."); } From 40397eb6266c4fe24406cff74440a1c9797d7c43 Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Thu, 4 May 2017 10:08:21 -0500 Subject: [PATCH 28/57] Remove rply from PlyWriter. --- io/PlyWriter.cpp | 258 ++++++++++++++++++++++-------------------- io/PlyWriter.hpp | 63 ++++++++--- pdal/PDALUtils.hpp | 4 +- pdal/Stage.cpp | 16 ++- pdal/Stage.hpp | 18 +++ pdal/util/OStream.hpp | 88 +++++++++++++- 6 files changed, 307 insertions(+), 140 deletions(-) diff --git a/io/PlyWriter.cpp b/io/PlyWriter.cpp index 278aae28ee..936ac19740 100644 --- a/io/PlyWriter.cpp +++ b/io/PlyWriter.cpp @@ -38,38 +38,11 @@ #include #include +#include #include namespace pdal { -namespace -{ - -void createErrorCallback(p_ply ply, const char* message) -{ - throw PlyWriter::error(message); -} - - -e_ply_type getPlyType(Dimension::Type type) -{ - static std::map types = - { - { Dimension::Type::Unsigned8, PLY_UINT8 }, - { Dimension::Type::Signed8, PLY_INT8 }, - { Dimension::Type::Unsigned16, PLY_UINT16 }, - { Dimension::Type::Signed16, PLY_INT16 }, - { Dimension::Type::Unsigned32, PLY_UINT32 }, - { Dimension::Type::Signed32, PLY_INT32 }, - { Dimension::Type::Float, PLY_FLOAT32 }, - { Dimension::Type::Double, PLY_FLOAT64 } - }; - - return types[type]; -} - -} // unnamed namespace - static PluginInfo const s_info = PluginInfo( "writers.ply", @@ -83,51 +56,18 @@ std::string PlyWriter::getName() const { return s_info.name; } PlyWriter::PlyWriter() - : m_ply(nullptr) - , m_pointCollector(nullptr) - , m_storageMode(PLY_DEFAULT) {} void PlyWriter::addArgs(ProgramArgs& args) { args.add("filename", "Output filename", m_filename).setPositional(); - args.add("storage_mode", "PLY Storage mode", m_storageModeSpec, "default"); + 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); } -void PlyWriter::initialize() -{ - std::string storageMode(m_storageModeSpec); - storageMode = Utils::tolower(storageMode); - - if (storageMode == "ascii") - { - m_storageMode = PLY_ASCII; - } - else if (storageMode == "little endian") - { - m_storageMode = PLY_LITTLE_ENDIAN; - } - else if (storageMode == "big endian") - { - m_storageMode = PLY_BIG_ENDIAN; - } - else if (storageMode == "default") - { - m_storageMode = PLY_DEFAULT; - } - else - { - throwError("Unknown storage mode '" + m_storageModeSpec + - "'. Known storage modes are: 'ascii', 'little endian', " - "'big endian', and 'default'"); - } -} - - void PlyWriter::prepared(PointTableRef table) { if (m_dimNames.size()) @@ -150,103 +90,175 @@ void PlyWriter::prepared(PointTableRef table) } -void PlyWriter::ready(PointTableRef table) +std::string PlyWriter::getType(Dimension::Type type) const { + static std::map types = + { + { Dimension::Type::Signed8, "int8" }, + { Dimension::Type::Unsigned8, "uint8" }, + { Dimension::Type::Signed16, "int16" }, + { Dimension::Type::Unsigned16, "uint16" }, + { Dimension::Type::Signed32, "int32" }, + { Dimension::Type::Unsigned32, "uint32" }, + { Dimension::Type::Float, "float32" }, + { Dimension::Type::Double, "float64" } + }; + try { - m_ply = ply_create(m_filename.c_str(), m_storageMode, - createErrorCallback, 0, nullptr); - if (!m_ply) - throwError("Could not open file '" + m_filename + "for writing."); + return types.at(type); } - catch(const error& err) + catch (std::out_of_range) { - throwError(err.what()); + throwError("Can't write dimension of type '" + + Dimension::interpretationName(type) + "'."); } - m_pointCollector.reset(new PointView(table)); + return ""; } -void PlyWriter::write(const PointViewPtr data) +void PlyWriter::writeHeader(PointLayoutPtr layout) const { - if (m_faces && m_pointCollector->size()) - throwError("Can't output faces for pipeline with multiple " - "point views."); - m_pointCollector->append(*data); + *m_stream << "ply" << std::endl; + *m_stream << "format " << m_format << " 1.0" << std::endl; + *m_stream << "comment Generated by PDAL" << std::endl; + *m_stream << "element vertex " << pointCount() << std::endl; + + auto ni = m_dimNames.begin(); + for (auto dim : m_dims) + { + std::string name = *ni++; + std::string typeString = getType(layout->dimType(dim)); + *m_stream << "property " << typeString << " " << name << std::endl; + } if (m_faces) { - m_mesh = data->mesh(); - if (!m_mesh) - throwError("Can't write mesh faces. No mesh has been generated."); + *m_stream << "element face " << faceCount() << std::endl; + *m_stream << "property list uint8 uint32 vertex_indices" << std::endl; } - if (data->size() > std::numeric_limits::max()) + *m_stream << "end_header" << std::endl; +} + + +void PlyWriter::ready(PointTableRef table) +{ + if (pointCount() > std::numeric_limits::max()) throwError("Can't write PLY file. Only " + std::to_string(std::numeric_limits::max()) + " points supported."); + + m_stream = Utils::createFile(m_filename, true); + writeHeader(table.layout()); } -void PlyWriter::done(PointTableRef table) +void PlyWriter::write(const PointViewPtr data) { - try - { - if (!ply_add_element(m_ply, "vertex", m_pointCollector->size())) - throwError("Could not add vertex element"); + m_views.push_back(data); +} - auto ni = m_dimNames.begin(); - for (auto dim : m_dims) { - std::string name = *ni++; - e_ply_type plyType = getPlyType(table.layout()->dimType(dim)); - if (!ply_add_scalar_property(m_ply, name.c_str(), plyType)) - throwError("Could not add scalar property '" + name + "'"); + +void PlyWriter::writeValue(PointRef& point, Dimension::Id dim, + Dimension::Type type) +{ + if (m_format == Format::ASCII) + { + double d = point.getFieldAs(dim); + *m_stream << d; } - if (m_faces) + else if (m_format == Format::BINARY_LE) + { + OLeStream out(m_stream); + Everything e; + point.getField((char *)&e, dim, type); + Utils::insertDim(out, type, e); + } + else if (m_format == Format::BINARY_BE) + { + OBeStream out(m_stream); + Everything e; + point.getField((char *)&e, dim, type); + Utils::insertDim(out, type, e); + } +} + + +void PlyWriter::writePoint(PointRef& point, PointLayoutPtr layout) +{ + for (auto it = m_dims.begin(); it != m_dims.end();) { - ply_add_element(m_ply, "face", m_mesh->size()); - ply_add_list_property(m_ply, "vertex_indices", PLY_UINT8, PLY_UINT32); + Dimension::Id dim = *it; + writeValue(point, dim, layout->dimType(dim)); + ++it; + if (m_format == Format::ASCII && it != m_dims.end()) + *m_stream << " "; } + if (m_format == Format::ASCII) + *m_stream << std::endl; +} - if (!ply_add_comment(m_ply, "Generated by PDAL")) - throwError("Could not add comment"); - if (!ply_write_header(m_ply)) - throwError("Could not write ply header"); +void PlyWriter::writeTriangle(const Triangle& t, size_t offset) +{ + if (m_format == Format::ASCII) + { + *m_stream << "3 " << (t.m_a + offset) << " " << + (t.m_b + offset) << " " << (t.m_c + offset) << std::endl; + } + else if (m_format == Format::BINARY_LE) + { + OLeStream out(m_stream); + unsigned char count = 3; + uint32_t a = (uint32_t)(t.m_a + offset); + uint32_t b = (uint32_t)(t.m_b + offset); + uint32_t c = (uint32_t)(t.m_c + offset); + out << count << a << b << c; + } + else if (m_format == Format::BINARY_BE) + { + OBeStream out(m_stream); + unsigned char count = 3; + uint32_t a = (uint32_t)(t.m_a + offset); + uint32_t b = (uint32_t)(t.m_b + offset); + uint32_t c = (uint32_t)(t.m_c + offset); + out << count << a << b << c; + } +} + - for (PointId index = 0; index < m_pointCollector->size(); ++index) +// Deferring write until this time allows both points and faces from multiple +// point views to be written. +void PlyWriter::done(PointTableRef table) +{ + for (auto& v : m_views) { - auto ni = m_dimNames.begin(); - auto di = m_dims.begin(); - for (; di != m_dims.end(); ++di, ++ni) + PointRef point(*v, 0); + for (PointId idx = 0; idx < v->size(); ++idx) { - double value = m_pointCollector->getFieldAs(*di, index); - if (!ply_write(m_ply, value)) - throwError("Error writing dimension '" + *ni + "' of point " - "number " + Utils::toString(index) + "."); + point.setPointId(idx); + writePoint(point, table.layout()); } } if (m_faces) { - for (PointId id = 0; id < m_mesh->size(); id++) + PointId offset = 0; + for (auto& v : m_views) { - const Triangle& t = (*m_mesh)[id]; - if (!(ply_write(m_ply, 3) && - ply_write(m_ply, t.m_a) && - ply_write(m_ply, t.m_b) && - ply_write(m_ply, t.m_c))) - throwError("Error writing face number " + - Utils::toString(id) + "."); + TriangularMesh *mesh = v->mesh(); + if (mesh) + { + for (size_t id = 0; id < mesh->size(); ++id) + { + const Triangle& t = (*mesh)[id]; + writeTriangle(t, offset); + } + } + offset += v->size(); } } - - if (!ply_close(m_ply)) - throwError("Error closing file"); - } - catch (const error& err) - { - throwError(err.what()); - } - + Utils::closeFile(m_stream); + m_stream = nullptr; getMetadata().addList("filename", m_filename); } -} +} // namespace pdal diff --git a/io/PlyWriter.hpp b/io/PlyWriter.hpp index f8b8fc6a30..02941cd737 100644 --- a/io/PlyWriter.hpp +++ b/io/PlyWriter.hpp @@ -32,8 +32,6 @@ * OF SUCH DAMAGE. ****************************************************************************/ -#include - #include #include #include @@ -43,18 +41,19 @@ extern "C" PF_ExitFunc PlyWriter_InitPlugin(); namespace pdal { -class TriangularMesh; + +class Triangle; class PDAL_DLL PlyWriter : public Writer { public: - struct error : public std::runtime_error + enum class Format { - error(const std::string& e) : std::runtime_error(e) - {} + ASCII, + BINARY_LE, + BINARY_BE }; -public: static void * create(); static int32_t destroy(void *); std::string getName() const; @@ -63,21 +62,59 @@ class PDAL_DLL PlyWriter : public Writer private: virtual void addArgs(ProgramArgs& args); - virtual void initialize(); virtual void prepared(PointTableRef table); virtual void ready(PointTableRef table); virtual void write(const PointViewPtr data); virtual void done(PointTableRef table); + std::string getType(Dimension::Type type) const; + void writeHeader(PointLayoutPtr layout) const; + void writeValue(PointRef& point, Dimension::Id dim, Dimension::Type type); + void writePoint(PointRef& point, PointLayoutPtr layout); + void writeTriangle(const Triangle& t, size_t offset); + + std::ostream *m_stream; std::string m_filename; - p_ply m_ply; - PointViewPtr m_pointCollector; - std::string m_storageModeSpec; - e_ply_storage_mode m_storageMode; + Format m_format; bool m_faces; StringList m_dimNames; Dimension::IdList m_dims; - TriangularMesh *m_mesh; + std::vector m_views; }; +inline std::istream& operator>>(std::istream& in, PlyWriter::Format& f) +{ + std::string s; + std::getline(in, s); + Utils::trim(s); + Utils::tolower(s); + if (s == "ascii" || s == "default") + f = PlyWriter::Format::ASCII; + else if (s == "little endian" || s == "binary_little_endian") + f = PlyWriter::Format::BINARY_LE; + else if (s == "big endian" || s == "binary_big_endian") + f = PlyWriter::Format::BINARY_BE; + else + in.setstate(std::ios_base::failbit); + return in; +} + + +inline std::ostream& operator<<(std::ostream& out, const PlyWriter::Format& f) +{ + switch (f) + { + case PlyWriter::Format::ASCII: + out << "ascii"; + break; + case PlyWriter::Format::BINARY_LE: + out << "binary_little_endian"; + break; + case PlyWriter::Format::BINARY_BE: + out << "binary_big_endian"; + break; + } + return out; +} + } diff --git a/pdal/PDALUtils.hpp b/pdal/PDALUtils.hpp index 90af65d1bd..7fd7749ebc 100644 --- a/pdal/PDALUtils.hpp +++ b/pdal/PDALUtils.hpp @@ -155,8 +155,8 @@ inline Everything extractDim(IN& ext, Dimension::Type type) } -inline void insertDim(Inserter& ins, Dimension::Type type, - const Everything& e) +template +inline void insertDim(OUT& ins, Dimension::Type type, const Everything& e) { using Type = Dimension::Type; diff --git a/pdal/Stage.cpp b/pdal/Stage.cpp index 76971a483c..d6ccc06ea0 100644 --- a/pdal/Stage.cpp +++ b/pdal/Stage.cpp @@ -49,7 +49,8 @@ namespace pdal { -Stage::Stage() : m_progressFd(-1), m_debug(false), m_verbose(0) +Stage::Stage() : m_progressFd(-1), m_debug(false), m_verbose(0), + m_pointCount(0), m_faceCount(0) {} @@ -178,6 +179,17 @@ PointViewSet Stage::execute(PointTableRef table) table.addSpatialReference((*it)->spatialReference()); gdal::ErrorHandler::getGlobalErrorHandler().set(m_log, m_debug); + // Count the number of views and the number of points and faces so they're + // available to stages. + m_pointCount = 0; + m_faceCount = 0; + for (auto const& it : views) + { + m_pointCount += it->size(); + auto m = it->mesh(); + if (m) + m_faceCount += m->size(); + } // Do the ready operation and then start running all the views // through the stage. ready(table); @@ -205,6 +217,8 @@ PointViewSet Stage::execute(PointTableRef table) } l_done(table); popLogLeader(); + m_pointCount = 0; + m_faceCount = 0; return outViews; } diff --git a/pdal/Stage.hpp b/pdal/Stage.hpp index d3d5c3ab03..55fd2bd2f6 100644 --- a/pdal/Stage.hpp +++ b/pdal/Stage.hpp @@ -321,6 +321,22 @@ class PDAL_DLL Stage void setSpatialReference(MetadataNode& m, SpatialReference const&); void addSpatialReferenceArg(ProgramArgs& args); void throwError(const std::string& s) const; + /** + Return the point count of all point views at the start of execution. + Only valid during execute(). + + \return Total number of points in all point views being executed. + */ + point_count_t pointCount() const + { return m_pointCount; } + /** + Return the count of faces in all primary meshes for all point views. + Only valid during execute(). + + \return Total number of faces in all point views being executed. + */ + point_count_t faceCount() const + { return m_faceCount; } private: bool m_debug; @@ -336,6 +352,8 @@ class PDAL_DLL Stage // bind the user_data argument that is essentially a comment in pipeline // files. std::string m_userDataJSON; + point_count_t m_pointCount; + point_count_t m_faceCount; Stage& operator=(const Stage&); // not implemented Stage(const Stage&); // not implemented diff --git a/pdal/util/OStream.hpp b/pdal/util/OStream.hpp index 9b2fc3b2c6..ca54caf662 100644 --- a/pdal/util/OStream.hpp +++ b/pdal/util/OStream.hpp @@ -206,6 +206,92 @@ class OLeStream : public OStream } }; + +/// Stream wrapper for output of binary data that converts from host ordering +/// to big endian format +class OBeStream : public OStream +{ +public: + PDAL_DLL OBeStream() + {} + PDAL_DLL OBeStream(const std::string& filename) : OStream(filename) + {} + PDAL_DLL OBeStream(std::ostream *stream) : OStream(stream) + {} + + PDAL_DLL OBeStream& operator << (uint8_t v) + { + m_stream->put((char)v); + return *this; + } + + PDAL_DLL OBeStream& operator << (int8_t v) + { + m_stream->put((char)v); + return *this; + } + + PDAL_DLL OBeStream& operator << (uint16_t v) + { + v = htobe16(v); + m_stream->write((char *)&v, sizeof(v)); + return *this; + } + + PDAL_DLL OBeStream& operator << (int16_t v) + { + v = (int16_t)htobe16((uint16_t)v); + m_stream->write((char *)&v, sizeof(v)); + return *this; + } + + PDAL_DLL OBeStream& operator << (uint32_t v) + { + v = htobe32(v); + m_stream->write((char *)&v, sizeof(v)); + return *this; + } + + PDAL_DLL OBeStream& operator << (int32_t v) + { + v = (int32_t)htobe32((uint32_t)v); + m_stream->write((char *)&v, sizeof(v)); + return *this; + } + + PDAL_DLL OBeStream& operator << (uint64_t v) + { + v = htobe64(v); + m_stream->write((char *)&v, sizeof(v)); + return *this; + } + + PDAL_DLL OBeStream& operator << (int64_t v) + { + v = (int64_t)htobe64((uint64_t)v); + m_stream->write((char *)&v, sizeof(v)); + return *this; + } + + PDAL_DLL OBeStream& operator << (float v) + { + uint32_t tmp(0); + std::memcpy(&tmp, &v, sizeof(v)); + tmp = htobe32(tmp); + m_stream->write((char *)&tmp, sizeof(tmp)); + return *this; + } + + PDAL_DLL OBeStream& operator << (double v) + { + uint64_t tmp(0); + std::memcpy(&tmp, &v, sizeof(v)); + tmp = htobe64(tmp); + m_stream->write((char *)&tmp, sizeof(tmp)); + return *this; + } +}; + /// Stream position marker with rewinding/reset support. class OStreamMarker { @@ -226,7 +312,7 @@ class OStreamMarker private: std::streampos m_pos; OStream& m_stream; - + OStreamMarker(const OStreamMarker&); // not implemented OStreamMarker& operator=(const OStreamMarker&); // not implemented }; From b3729504671ebea8d66b17bcd53b00cbb618aa9a Mon Sep 17 00:00:00 2001 From: Andrew Bell Date: Thu, 4 May 2017 10:29:28 -0500 Subject: [PATCH 29/57] Remove rply. Use came case for ply format. --- CMakeLists.txt | 1 - cmake/rply.cmake | 1 - io/PlyReader.cpp | 18 +- io/PlyReader.hpp | 6 +- io/PlyWriter.cpp | 18 +- io/PlyWriter.hpp | 18 +- vendor/rply/LICENSE | 20 - vendor/rply/etc/convert.c | 130 --- vendor/rply/etc/dump.c | 44 - vendor/rply/etc/input.ply | 16 - vendor/rply/etc/sconvert.c | 64 -- vendor/rply/manual/manual.html | 1137 --------------------- vendor/rply/manual/reference.css | 54 - vendor/rply/manual/rply.png | Bin 6232 -> 0 bytes vendor/rply/rply.c | 1616 ------------------------------ vendor/rply/rply.h | 378 ------- vendor/rply/rplyfile.h | 68 -- 17 files changed, 30 insertions(+), 3559 deletions(-) delete mode 100644 cmake/rply.cmake delete mode 100644 vendor/rply/LICENSE delete mode 100644 vendor/rply/etc/convert.c delete mode 100644 vendor/rply/etc/dump.c delete mode 100644 vendor/rply/etc/input.ply delete mode 100644 vendor/rply/etc/sconvert.c delete mode 100644 vendor/rply/manual/manual.html delete mode 100644 vendor/rply/manual/reference.css delete mode 100644 vendor/rply/manual/rply.png delete mode 100644 vendor/rply/rply.c delete mode 100644 vendor/rply/rply.h delete mode 100644 vendor/rply/rplyfile.h diff --git a/CMakeLists.txt b/CMakeLists.txt index a91cee141b..568316b2fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -119,7 +119,6 @@ include(${PDAL_CMAKE_DIR}/curl.cmake) if (CURL_FOUND) set(PDAL_ARBITER_LIB_NAME pdal_arbiter) endif() -include(${PDAL_CMAKE_DIR}/rply.cmake) include(${PDAL_CMAKE_DIR}/json.cmake) include(${PDAL_CMAKE_DIR}/libxml2.cmake) include(${PDAL_CMAKE_DIR}/dimension.cmake) diff --git a/cmake/rply.cmake b/cmake/rply.cmake deleted file mode 100644 index a242529055..0000000000 --- a/cmake/rply.cmake +++ /dev/null @@ -1 +0,0 @@ -set(RPLY_SRCS ${PDAL_VENDOR_DIR}/rply/rply.c) diff --git a/io/PlyReader.cpp b/io/PlyReader.cpp index 9df86df1e8..5e20100ab7 100644 --- a/io/PlyReader.cpp +++ b/io/PlyReader.cpp @@ -118,11 +118,11 @@ void PlyReader::extractFormat() word = nextWord(); if (word == "ascii") - m_format = Format::ASCII; + m_format = Format::Ascii; else if (word == "binary_big_endian") - m_format = Format::BINARY_BE; + m_format = Format::BinaryBe; else if (word == "binary_little_endian") - m_format = Format::BINARY_LE; + m_format = Format::BinaryLe; else throwError("Unrecognized PLY format: '" + word + "'."); @@ -325,19 +325,19 @@ bool PlyReader::readProperty(Property *prop, PointRef& point) void PlyReader::SimpleProperty::read(std::istream *stream, PlyReader::Format format, PointRef& point) { - if (format == Format::ASCII) + if (format == Format::Ascii) { double d; *stream >> d; point.setField(m_dim, d); } - else if (format == Format::BINARY_LE) + else if (format == Format::BinaryLe) { ILeStream in(stream); Everything e = Utils::extractDim(in, m_type); point.setField(m_dim, m_type, &e); } - else if (format == Format::BINARY_BE) + else if (format == Format::BinaryBe) { IBeStream in(stream); Everything e = Utils::extractDim(in, m_type); @@ -351,7 +351,7 @@ void PlyReader::SimpleProperty::read(std::istream *stream, void PlyReader::ListProperty::read(std::istream *stream, PlyReader::Format format, PointRef& point) { - if (format == Format::ASCII) + if (format == Format::Ascii) { size_t cnt; *stream >> cnt; @@ -360,7 +360,7 @@ void PlyReader::ListProperty::read(std::istream *stream, while (cnt--) *stream >> d; } - else if (format == Format::BINARY_LE) + else if (format == Format::BinaryLe) { ILeStream istream(stream); Everything e = Utils::extractDim(istream, m_countType); @@ -368,7 +368,7 @@ void PlyReader::ListProperty::read(std::istream *stream, cnt *= Dimension::size(m_listType); istream.seek(cnt, std::ios_base::cur); } - else if (format == Format::BINARY_BE) + else if (format == Format::BinaryBe) { IBeStream istream(stream); Everything e = Utils::extractDim(istream, m_countType); diff --git a/io/PlyReader.hpp b/io/PlyReader.hpp index 3b4d8657e1..9f5949caa9 100644 --- a/io/PlyReader.hpp +++ b/io/PlyReader.hpp @@ -62,9 +62,9 @@ class PDAL_DLL PlyReader : public Reader private: enum class Format { - ASCII, - BINARY_LE, - BINARY_BE + Ascii, + BinaryLe, + BinaryBe }; struct Property diff --git a/io/PlyWriter.cpp b/io/PlyWriter.cpp index 936ac19740..9bd9f045b8 100644 --- a/io/PlyWriter.cpp +++ b/io/PlyWriter.cpp @@ -62,7 +62,7 @@ PlyWriter::PlyWriter() void PlyWriter::addArgs(ProgramArgs& args) { args.add("filename", "Output filename", m_filename).setPositional(); - args.add("storage_mode", "PLY Storage Mode", m_format, Format::ASCII); + 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); } @@ -161,19 +161,19 @@ void PlyWriter::write(const PointViewPtr data) void PlyWriter::writeValue(PointRef& point, Dimension::Id dim, Dimension::Type type) { - if (m_format == Format::ASCII) + if (m_format == Format::Ascii) { double d = point.getFieldAs(dim); *m_stream << d; } - else if (m_format == Format::BINARY_LE) + else if (m_format == Format::BinaryLe) { OLeStream out(m_stream); Everything e; point.getField((char *)&e, dim, type); Utils::insertDim(out, type, e); } - else if (m_format == Format::BINARY_BE) + else if (m_format == Format::BinaryBe) { OBeStream out(m_stream); Everything e; @@ -190,22 +190,22 @@ void PlyWriter::writePoint(PointRef& point, PointLayoutPtr layout) Dimension::Id dim = *it; writeValue(point, dim, layout->dimType(dim)); ++it; - if (m_format == Format::ASCII && it != m_dims.end()) + if (m_format == Format::Ascii && it != m_dims.end()) *m_stream << " "; } - if (m_format == Format::ASCII) + if (m_format == Format::Ascii) *m_stream << std::endl; } void PlyWriter::writeTriangle(const Triangle& t, size_t offset) { - if (m_format == Format::ASCII) + if (m_format == Format::Ascii) { *m_stream << "3 " << (t.m_a + offset) << " " << (t.m_b + offset) << " " << (t.m_c + offset) << std::endl; } - else if (m_format == Format::BINARY_LE) + else if (m_format == Format::BinaryLe) { OLeStream out(m_stream); unsigned char count = 3; @@ -214,7 +214,7 @@ void PlyWriter::writeTriangle(const Triangle& t, size_t offset) uint32_t c = (uint32_t)(t.m_c + offset); out << count << a << b << c; } - else if (m_format == Format::BINARY_BE) + else if (m_format == Format::BinaryBe) { OBeStream out(m_stream); unsigned char count = 3; diff --git a/io/PlyWriter.hpp b/io/PlyWriter.hpp index 02941cd737..7b43e8885f 100644 --- a/io/PlyWriter.hpp +++ b/io/PlyWriter.hpp @@ -49,9 +49,9 @@ class PDAL_DLL PlyWriter : public Writer public: enum class Format { - ASCII, - BINARY_LE, - BINARY_BE + Ascii, + BinaryLe, + BinaryBe }; static void * create(); @@ -89,11 +89,11 @@ inline std::istream& operator>>(std::istream& in, PlyWriter::Format& f) Utils::trim(s); Utils::tolower(s); if (s == "ascii" || s == "default") - f = PlyWriter::Format::ASCII; + f = PlyWriter::Format::Ascii; else if (s == "little endian" || s == "binary_little_endian") - f = PlyWriter::Format::BINARY_LE; + f = PlyWriter::Format::BinaryLe; else if (s == "big endian" || s == "binary_big_endian") - f = PlyWriter::Format::BINARY_BE; + f = PlyWriter::Format::BinaryBe; else in.setstate(std::ios_base::failbit); return in; @@ -104,13 +104,13 @@ inline std::ostream& operator<<(std::ostream& out, const PlyWriter::Format& f) { switch (f) { - case PlyWriter::Format::ASCII: + case PlyWriter::Format::Ascii: out << "ascii"; break; - case PlyWriter::Format::BINARY_LE: + case PlyWriter::Format::BinaryLe: out << "binary_little_endian"; break; - case PlyWriter::Format::BINARY_BE: + case PlyWriter::Format::BinaryBe: out << "binary_big_endian"; break; } diff --git a/vendor/rply/LICENSE b/vendor/rply/LICENSE deleted file mode 100644 index ca893cb03d..0000000000 --- a/vendor/rply/LICENSE +++ /dev/null @@ -1,20 +0,0 @@ -RPly 1.1.4 license -Copyright © 2003-2015 Diego Nehab. - -Permission is hereby granted, free of charge, to any person obtaining a -copy of this software and associated documentation files (the "Software"), -to deal in the Software without restriction, including without limitation -the rights to use, copy, modify, merge, publish, distribute, sublicense, -and/or sell copies of the Software, and to permit persons to whom the -Software is furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in -all copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER -DEALINGS IN THE SOFTWARE. diff --git a/vendor/rply/etc/convert.c b/vendor/rply/etc/convert.c deleted file mode 100644 index e87ac09caf..0000000000 --- a/vendor/rply/etc/convert.c +++ /dev/null @@ -1,130 +0,0 @@ -#include -#include -#include -#include -#include "rply.h" - -/* internal function prototypes */ -static void error(const char *fmt, ...); -static void help(void); -static void parse_arguments(int argc, char **argv, - e_ply_storage_mode *storage_mode, - const char **iname, const char **oname); -static void setup_callbacks(p_ply iply, p_ply oply); - -/* given a format mode, an input file name and an output file name, - * convert input file to output in given format mode */ -int main(int argc, char **argv) { - const char *value = NULL; - e_ply_storage_mode storage_mode = PLY_LITTLE_ENDIAN; - const char *iname = NULL, *oname = NULL; - p_ply iply = NULL, oply = NULL; - /* parse command line arguments */ - parse_arguments(argc, argv, &storage_mode, &iname, &oname); - /* open input file and make sure we parsed its header */ - iply = ply_open(iname, NULL, 0, NULL); - if (!iply) error("Unable to open file '%s'", iname); - if (!ply_read_header(iply)) error("Failed reading '%s' header", iname); - /* create output file */ - oply = ply_create(oname, storage_mode, NULL, 0, NULL); - if (!oply) error("Unable to create file '%s'", oname); - /* create elements and properties in output file and - * setup callbacks for them in input file */ - setup_callbacks(iply, oply); - /* pass comments and obj_infos from input to output */ - value = NULL; - while ((value = ply_get_next_comment(iply, value))) - if (!ply_add_comment(oply, value)) - error("Failed adding comments"); - value = NULL; - while ((value = ply_get_next_obj_info(iply, value))) - if (!ply_add_obj_info(oply, value)) - error("Failed adding comments"); - /* write output header */ - if (!ply_write_header(oply)) error("Failed writing '%s' header", oname); - /* read input file generating callbacks that pass data to output file */ - if (!ply_read(iply)) error("Conversion failed"); - /* close up, we are done */ - if (!ply_close(iply)) error("Error closing file '%s'", iname); - if (!ply_close(oply)) error("Error closing file '%s'", oname); - return 0; -} - -/* prints an error message and exits */ -static void error(const char *fmt, ...) { - va_list ap; - va_start(ap, fmt); - vfprintf(stderr, fmt, ap); - va_end(ap); - fprintf(stderr, "\n"); - exit(1); -} - -/* prints the help message and exits */ -static void help(void) { - error("Usage:\n" - " convert