diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h index 7de04ff26a9..7d2bb4842bd 100644 --- a/BGL/include/CGAL/boost/graph/parameters_interface.h +++ b/BGL/include/CGAL/boost/graph/parameters_interface.h @@ -152,6 +152,8 @@ CGAL_add_named_parameter(inspector_t, inspector, inspector) CGAL_add_named_parameter(logger_t, logger, logger) CGAL_add_named_parameter(pointmatcher_config_t, pointmatcher_config, pointmatcher_config) CGAL_add_named_parameter(adjacencies_t, adjacencies, adjacencies) +CGAL_add_named_parameter(scan_angle_t, scan_angle_map, scan_angle_map) +CGAL_add_named_parameter(scanline_id_t, scanline_id_map, scanline_id_map) // List of named parameters used in Surface_mesh_approximation package CGAL_add_named_parameter(verbose_level_t, verbose_level, verbose_level) diff --git a/BGL/test/BGL/test_cgal_bgl_named_params.cpp b/BGL/test/BGL/test_cgal_bgl_named_params.cpp index 91ebf6224db..18db63159e4 100644 --- a/BGL/test/BGL/test_cgal_bgl_named_params.cpp +++ b/BGL/test/BGL/test_cgal_bgl_named_params.cpp @@ -268,6 +268,8 @@ void test(const NamedParameters& np) check_same_type<9032>(get_parameter(np, CGAL::internal_np::inspector)); check_same_type<9033>(get_parameter(np, CGAL::internal_np::logger)); check_same_type<9034>(get_parameter(np, CGAL::internal_np::maximum_normal_deviation)); + check_same_type<9035>(get_parameter(np, CGAL::internal_np::scan_angle_map)); + check_same_type<9036>(get_parameter(np, CGAL::internal_np::scanline_id_map)); } int main() @@ -396,6 +398,8 @@ int main() .inspector(A<9032>(9032)) .logger(A<9033>(9033)) .maximum_normal_deviation(A<9034>(9034)) + .scan_angle_map(A<9035>(9035)) + .scanline_id_map(A<9036>(9036)) .maximum_number_of_faces(A<78910>(78910)) ); return EXIT_SUCCESS; diff --git a/Documentation/doc/Documentation/Tutorials/Tutorial_reconstruction.txt b/Documentation/doc/Documentation/Tutorials/Tutorial_reconstruction.txt index fb60465bbc0..7a3986f3158 100644 --- a/Documentation/doc/Documentation/Tutorials/Tutorial_reconstruction.txt +++ b/Documentation/doc/Documentation/Tutorials/Tutorial_reconstruction.txt @@ -184,9 +184,21 @@ PCA is faster but jet is more accurate in the presence of high curvatures. These function only estimates the _direction_ of the normals, not their orientation (the orientation of the vectors might not be locally consistent). To properly orient the normals, the -function `mst_orient_normals()` can be used. Notice that it can also -be used directly on input normals if their orientation is not -consistent. +following functions can be used: + +- `mst_orient_normals()` +- `scanline_orient_normals()` + +The first one uses a _minimum spanning tree_ to consistently propagate +the orientation of normals in an increasingly large neighborhood. In +the case of data with many sharp features and occlusions (which are +common in airborne LIDAR data, for example), the second algorithm may +produce better results: it takes advantage of point clouds which are +ordered into scanlines to estimate the line of sight of each point and +thus to orient normals accordingly. + +Notice that these can also be used directly on input normals if their +orientation is not consistent. \snippet Poisson_surface_reconstruction_3/tutorial_example.cpp Normal estimation diff --git a/Installation/CHANGES.md b/Installation/CHANGES.md index cd1400b49c9..cd587ece395 100644 --- a/Installation/CHANGES.md +++ b/Installation/CHANGES.md @@ -54,6 +54,9 @@ Release date: December 2020 - Added an optional range parameter to `CGAL::Polygon_mesh_processing::stitch_borders()`, which can be used to specify which boundary cycles are eligible for stitching. +### [Point Set Processing](https://doc.cgal.org/5.2/Manual/packages.html#PkgPointSetProcessing3) +- Added a function [`CGAL::scanline_orient_normals()`](https://doc.cgal.org/5.2/Point_set_processing_3/group__PkgPointSetProcessing3Algorithms.html#ga221d4efde44f42aefe153cb927138efe) that orients a point cloud by estimating a line of sight. + ### [Classification](https://doc.cgal.org/5.2/Manual/packages.html#PkgClassification) - **Breaking change**: new IO format for `ETHZ::Random_Forest` (a diff --git a/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt b/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt index 3d482fdc4f2..32e0e9f61c2 100644 --- a/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt +++ b/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt @@ -66,6 +66,7 @@ format. - `CGAL::jet_estimate_normals()` - `CGAL::pca_estimate_normals()` - `CGAL::mst_orient_normals()` +- `CGAL::scanline_orient_normals()` - `CGAL::edge_aware_upsample_point_set()` - `CGAL::compute_vcm()` - `CGAL::vcm_estimate_normals()` @@ -101,4 +102,3 @@ format. - `CGAL::make_las_point_writer()` */ - diff --git a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt index 514afa2b06c..9d7f9c48d89 100644 --- a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt +++ b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt @@ -832,6 +832,8 @@ of nearest neighbors or a fixed spherical neighborhood radius. \section Point_set_processing_3NormalOrientation Normal Orientation +\subsection Point_set_processing_3Mst_orient_normals Minimum Spanning Tree + Function `mst_orient_normals()` orients the normals of a set of points with unoriented normals using the method described by Hoppe et al. in Surface reconstruction from unorganized points \cgalCite{cgal:hddms-srup-92}. @@ -846,7 +848,7 @@ the normals which cannot be successfully oriented. Normal orientation of a sampled cube surface. Left: unoriented normals. Right: orientation of right face normals is propagated to bottom face. \cgalFigureEnd -\subsection Point_set_processing_3Example_normals Example +\subsubsection Point_set_processing_3Example_normals Example The following example reads a point set from a file, estimates the normals through PCA (either over the 18 nearest neighbors or using a @@ -854,8 +856,31 @@ spherical neighborhood radius of twice the average spacing) and orients the normals: \cgalExample{Point_set_processing_3/normals_example.cpp} +\subsection Point_set_processing_3Scanline_orient_normals Scanline +The minimum spanning tree results can give suboptimal results on point +clouds with many sharp features and occlusions, which typically +happens on airborne acquired urban datasets. +`scanline_orient_normals()` is an alternative method specialized for +point sets which are ordered in scanline aligned on the XY-plane. It +can take advantage of LAS properties provided by some LIDAR scanner +and is the best choice of normal orientation when dealing with 2.5D +urban scenes. + +\cgalFigureBegin{Point_set_processing_3figmst_scanline_orient_normals,scanline_orient_normals.png} +Normal orientation of a LIDAR scanline. The point cloud is a typical +airborne LIDAR input, sampling a building without normal information +and with many occlusions (especially on vertical walls). +\cgalFigureEnd + +\subsubsection Point_set_processing_3Example_scanline_normals Example + +The following example reads a point set from a LAS file, estimates the +normals through Jet Fitting and outputs in PLY format the orientation +results of all the variants of `scanline_orient_normals()`: + +\cgalExample{Point_set_processing_3/orient_scanlines_example.cpp} \section Point_set_processing_3Upsampling Upsampling diff --git a/Point_set_processing_3/doc/Point_set_processing_3/examples.txt b/Point_set_processing_3/doc/Point_set_processing_3/examples.txt index 2374e1ea02f..ea670af1b62 100644 --- a/Point_set_processing_3/doc/Point_set_processing_3/examples.txt +++ b/Point_set_processing_3/doc/Point_set_processing_3/examples.txt @@ -16,6 +16,7 @@ \example Point_set_processing_3/hierarchy_simplification_example.cpp \example Point_set_processing_3/jet_smoothing_example.cpp \example Point_set_processing_3/normals_example.cpp +\example Point_set_processing_3/orient_scanlines_example.cpp \example Point_set_processing_3/wlop_simplify_and_regularize_point_set_example.cpp \example Point_set_processing_3/bilateral_smooth_point_set_example.cpp \example Point_set_processing_3/edge_aware_upsample_point_set_example.cpp diff --git a/Point_set_processing_3/doc/Point_set_processing_3/fig/scanline_orient_normals.png b/Point_set_processing_3/doc/Point_set_processing_3/fig/scanline_orient_normals.png new file mode 100644 index 00000000000..2da96f7ea3f Binary files /dev/null and b/Point_set_processing_3/doc/Point_set_processing_3/fig/scanline_orient_normals.png differ diff --git a/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt b/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt index 83253368016..1e55a2dda63 100644 --- a/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt +++ b/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt @@ -59,18 +59,24 @@ if ( CGAL_FOUND ) target_link_libraries(${target} ${CGAL_libs}) endforeach() + # Use Eigen + find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater) + include(CGAL_Eigen_support) + find_package(LASLIB) include(CGAL_LASLIB_support) if (TARGET CGAL::LASLIB_support) add_executable( read_las_example "read_las_example.cpp" ) target_link_libraries(read_las_example ${CGAL_libs} CGAL::LASLIB_support) + if (TARGET CGAL::Eigen_support) + add_executable( orient_scanlines_example "orient_scanlines_example.cpp" ) + target_link_libraries(orient_scanlines_example ${CGAL_libs} + CGAL::Eigen_support CGAL::LASLIB_support) + endif() else() message(STATUS "NOTICE : the LAS reader test requires LASlib and will not be compiled.") endif() - # Use Eigen - find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater) - include(CGAL_Eigen_support) if (TARGET CGAL::Eigen_support) set(CGAL_libs ${CGAL_libs} CGAL::Eigen_support) diff --git a/Point_set_processing_3/examples/Point_set_processing_3/data/urban.las b/Point_set_processing_3/examples/Point_set_processing_3/data/urban.las new file mode 100644 index 00000000000..57da02b6301 Binary files /dev/null and b/Point_set_processing_3/examples/Point_set_processing_3/data/urban.las differ diff --git a/Point_set_processing_3/examples/Point_set_processing_3/orient_scanlines_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/orient_scanlines_example.cpp new file mode 100644 index 00000000000..90bb4a11917 --- /dev/null +++ b/Point_set_processing_3/examples/Point_set_processing_3/orient_scanlines_example.cpp @@ -0,0 +1,88 @@ +#include +#include +#include +#include +#include + +using Kernel = CGAL::Simple_cartesian; +using Point_3 = Kernel::Point_3; +using Vector_3 = Kernel::Vector_3; +using Point_with_info = std::tuple; +using Point_map = CGAL::Nth_of_tuple_property_map<0, Point_with_info>; +using Normal_map = CGAL::Nth_of_tuple_property_map<1, Point_with_info>; +using Scan_angle_map = CGAL::Nth_of_tuple_property_map<2, Point_with_info>; +using Scanline_id_map = CGAL::Nth_of_tuple_property_map<3, Point_with_info>; + +void dump (const char* filename, const std::vector& points) +{ + std::ofstream ofile (filename, std::ios::binary); + CGAL::set_binary_mode(ofile); + CGAL::write_ply_points + (ofile, points, + CGAL::parameters::point_map (Point_map()). + normal_map (Normal_map())); + +} + +int main (int argc, char** argv) +{ + std::string fname (argc > 1 ? argv[1] : "data/urban.las"); + + std::vector points; + + std::cerr << "Reading input file " << fname << std::endl; + std::ifstream ifile (fname, std::ios::binary); + if (!ifile || + !CGAL::read_las_points_with_properties + (ifile, std::back_inserter (points), + CGAL::make_las_point_reader (Point_map()), + std::make_pair (Scan_angle_map(), + CGAL::LAS_property::Scan_angle()), + std::make_pair (Scanline_id_map(), + CGAL::LAS_property::Scan_direction_flag()))) + { + std::cerr << "Can't read " << fname << std::endl; + return EXIT_FAILURE; + } + + + std::cerr << "Estimating normals" << std::endl; + CGAL::jet_estimate_normals + (points, 12, + CGAL::parameters::point_map (Point_map()). + normal_map (Normal_map())); + + std::cerr << "Orienting normals using scan angle and direction flag" << std::endl; + CGAL::scanline_orient_normals + (points, + CGAL::parameters::point_map (Point_map()). + normal_map (Normal_map()). + scan_angle_map (Scan_angle_map()). + scanline_id_map (Scanline_id_map())); + dump("out_angle_and_flag.ply", points); + + std::cerr << "Orienting normals using scan direction flag only" << std::endl; + CGAL::scanline_orient_normals + (points, + CGAL::parameters::point_map (Point_map()). + normal_map (Normal_map()). + scanline_id_map (Scanline_id_map())); + dump("out_flag.ply", points); + + std::cerr << "Orienting normals using scan angle only" << std::endl; + CGAL::scanline_orient_normals + (points, + CGAL::parameters::point_map (Point_map()). + normal_map (Normal_map()). + scan_angle_map (Scan_angle_map())); + dump("out_angle.ply", points); + + std::cerr << "Orienting normals using no additional info" << std::endl; + CGAL::scanline_orient_normals + (points, + CGAL::parameters::point_map (Point_map()). + normal_map (Normal_map())); + dump("out_nothing.ply", points); + + return EXIT_SUCCESS; +} diff --git a/Point_set_processing_3/include/CGAL/scanline_orient_normals.h b/Point_set_processing_3/include/CGAL/scanline_orient_normals.h new file mode 100644 index 00000000000..41cfbc15bfa --- /dev/null +++ b/Point_set_processing_3/include/CGAL/scanline_orient_normals.h @@ -0,0 +1,564 @@ +// Copyright (c) 2020 GeometryFactory Sarl (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s) : Simon Giraudot + +#ifndef CGAL_SCANLINE_ORIENT_NORMALS_H +#define CGAL_SCANLINE_ORIENT_NORMALS_H + +#include + +#include + +#include +#include + +#include + +#include + +#if defined(CGAL_EIGEN3_ENABLED) || defined(DOXYGEN_RUNNING) + +#include +#include + +//#define CGAL_SCANLINE_ORIENT_VERBOSE + +namespace CGAL +{ + +/// \cond SKIP_IN_MANUAL +namespace Point_set_processing_3 +{ + +namespace internal +{ + +template +const Vector_3& vertical_vector() +{ + static Vector_3 v(0, 0, 1); + return v; +} + +template +bool is_end_of_scanline (Iterator scanline_begin, Iterator it, + PointMap, + ScanlineIDMap scanline_id_map, + const Tag_false&) // no fallback +{ + return (get (scanline_id_map, *scanline_begin) + != get (scanline_id_map, *it)); +} + +template +bool is_end_of_scanline (Iterator scanline_begin, Iterator it, + PointMap point_map, + ScanlineIDMap, + const Tag_true&) // fallback +{ + using Point_3 = typename boost::property_traits::value_type; + using Vector_3 = typename Kernel_traits::Kernel::Vector_3; + + if (std::distance (scanline_begin, it) < 3) + return false; + + Iterator n_minus_1 = it; -- n_minus_1; + Iterator n_minus_2 = n_minus_1; -- n_minus_2; + Iterator n_minus_3 = n_minus_2; -- n_minus_3; + + const Point_3& p_minus_1 = get(point_map, *n_minus_1); + const Point_3& p_minus_2 = get(point_map, *n_minus_2); + const Point_3& p_minus_3 = get(point_map, *n_minus_3); + + // End of scanline reached if inversion of direction + Vector_3 v32 (p_minus_3, p_minus_2); + v32 = Vector_3 (v32.x(), v32.y(), 0); + Vector_3 v21 (p_minus_2, p_minus_1); + v21 = Vector_3 (v21.x(), v21.y(), 0); + + return (v32 * v21 < 0); +} + +template +std::pair::value_type>::Kernel::Vector_3, + typename Kernel_traits::value_type>::Kernel::Vector_3> +scanline_base (Iterator begin, Iterator end, + PointMap point_map) +{ + using Point_3 = typename boost::property_traits::value_type; + using Kernel = typename Kernel_traits::Kernel; + using Vector_3 = typename Kernel::Vector_3; + + using Line_3 = typename Kernel::Line_3; + using Plane_3 = typename Kernel::Plane_3; + + const double limit = CGAL_PI * 30. / 180.; + + Line_3 pca2; + linear_least_squares_fitting_3 + (boost::make_transform_iterator + (begin, Property_map_to_unary_function(point_map)), + boost::make_transform_iterator + (end, Property_map_to_unary_function(point_map)), + pca2, Dimension_tag<0>()); + + Vector_3 pca_direction = pca2.to_vector(); + pca_direction = Vector_3 (pca_direction.x(), pca_direction.y(), 0); + pca_direction = pca_direction / CGAL::approximate_sqrt(pca_direction.squared_length()); + + Plane_3 pca3; + linear_least_squares_fitting_3 + (boost::make_transform_iterator + (begin, Property_map_to_unary_function(point_map)), + boost::make_transform_iterator + (end, Property_map_to_unary_function(point_map)), + pca3, Dimension_tag<0>()); + + Vector_3 orthogonal = pca3.orthogonal_vector(); + Vector_3 vertical = CGAL::cross_product (pca_direction, orthogonal); + if (vertical * vertical_vector() < 0) + vertical = -vertical; + + vertical = vertical / CGAL::approximate_sqrt(vertical.squared_length()); + + if (std::acos(vertical * vertical_vector()) < limit) + return std::make_pair (pca_direction, vertical); + + // if plane diverges from the vertical more than 30 degrees, then + // fallback to 0 0 1 vertical vector + + // Dummy begin->end vector version + Iterator last = end; -- last; + const Point_3& pbegin = get (point_map, *begin); + const Point_3& plast = get (point_map, *last); + Vector_3 direction (Point_3 (pbegin.x(), pbegin.y(), 0), + Point_3 (plast.x(), plast.y(), 0)); + direction = direction / CGAL::approximate_sqrt (direction.squared_length()); + + return std::make_pair (direction, vertical_vector()); +} + +template +bool normal_along_scanline_is_inverted +(const Vector_3& normal, const Vector_3& line_of_sight) +{ + return (line_of_sight * normal < 0); +} + +template +std::pair::value_type, bool> +estimate_scan_position (Iterator begin, Iterator end, PointMap point_map, + const std::vector& lines_of_sight) +{ + using Point_3 = typename boost::property_traits::value_type; + typedef Eigen::Matrix3d Matrix; + typedef Eigen::Vector3d Vector; + typedef Eigen::ConjugateGradient Solver; + + Matrix R; + R << 0, 0, 0, 0, 0, 0, 0, 0, 0; + Vector Q; + Q << 0, 0, 0; + + std::size_t idx = 0; + for (Iterator it = begin; it != end; ++ it, ++ idx) + { + const Point_3& p = get (point_map, *it); + const Vector_3& v = lines_of_sight[idx]; + + Vector n; + n << v.x(), v.y(), v.z(); + + Matrix I_nnt = Matrix::Identity() - n * n.transpose(); + Vector a; + a << p.x(), p.y(), p.z(); + + R += I_nnt; + Q += I_nnt * a; + } + + Solver solver(R); + + if (solver.info() != Eigen::Success) + return std::make_pair (ORIGIN, false); + + Vector p = solver.solve(Q); + if (solver.info() != Eigen::Success) + return std::make_pair (ORIGIN, false); + + return std::make_pair (Point_3(p(0), p(1), p(2)), true); +} + + +template +void orient_scanline (Iterator begin, Iterator end, + PointMap point_map, + NormalMap normal_map, + ScanAngleMap scan_angle_map, + const Tag_false&) // no fallback scan angle +{ + using Point_3 = typename boost::property_traits::value_type; + using Vector_3 = typename boost::property_traits::value_type; + + Vector_3 direction; + Vector_3 vertical; + std::tie (direction, vertical) + = scanline_base (begin, end, point_map); + + std::vector lines_of_sight; + lines_of_sight.reserve (std::distance (begin, end)); + + double mean_z = 0; + for (Iterator it = begin; it != end; ++ it) + { + double angle = CGAL_PI * static_cast(get (scan_angle_map, *it)) / 180.; + Vector_3 los = direction * std::sin(angle) + vertical * std::cos(angle); + lines_of_sight.push_back (los); + mean_z += get (point_map, *it).z(); + } + mean_z /= std::distance (begin, end); + +#ifdef CGAL_SCANORIENT_DUMP_RANDOM_SCANLINES + if (rand() % 1000 == 0 && std::distance(begin, end) > 10) + { + std::ofstream ofile ("scanline.polylines.txt"); + ofile.precision(18); + ofile << std::distance (begin, end); + for (Iterator it = begin; it != end; ++ it) + ofile << " " << get(point_map, *it); + ofile << std::endl; + + std::ofstream ofile2 ("base.polylines.txt"); + Point_3 orig = get (point_map, *(begin + std::distance(begin, end) / 2)); + double dist = CGAL::approximate_sqrt (CGAL::squared_distance + (get(point_map, *begin), orig)); + + ofile2.precision(18); + ofile2 << "2 " << orig << " " << orig + direction * dist << std::endl + << "2 " << orig << " " << orig + vertical * dist << std::endl; + } +#endif + + Point_3 scan_position; + bool solver_success; + std::tie (scan_position, solver_success) + = estimate_scan_position (begin, end, point_map, lines_of_sight); + + // If solver failed OR if scan position is detected under the + // scanline (obviously wrong) + if (!solver_success || scan_position.z() < mean_z) + { +#ifdef CGAL_SCANLINE_ORIENT_VERBOSE + if (!solver_success) + std::cerr << "Inverting because olver failed: "; + else + std::cerr << "Inverting because scanner under scanline: "; +#endif + + direction = -direction; + std::size_t idx = 0; + for (Iterator it = begin; it != end; ++ it, ++ idx) + { + double angle = CGAL_PI * static_cast(get (scan_angle_map, *it)) / 180.; + Vector_3 los = direction * std::sin(angle) + vertical * std::cos(angle); + lines_of_sight[idx] = los; + } + + std::tie (scan_position, solver_success) + = estimate_scan_position (begin, end, point_map, lines_of_sight); + +#ifdef CGAL_SCANLINE_ORIENT_VERBOSE + if (solver_success && scan_position.z() > mean_z) + std::cerr << "SOLVED" << std::endl; + else if (!solver_success) + std::cerr << "FAILED, solver failure" << std::endl; + else + std::cerr << "FAILED, scanner under scanline" << std::endl; +#endif + } + + std::size_t idx = 0; + for (Iterator it = begin; it != end; ++ it, ++ idx) + { + const Vector_3 los = lines_of_sight[idx]; + const Vector_3& normal = get (normal_map, *it); + if (normal_along_scanline_is_inverted (normal, los)) + put (normal_map, *it, -normal); + } +} + +template +void orient_scanline (Iterator begin, Iterator end, + PointMap point_map, + NormalMap normal_map, + ScanAngleMap, + const Tag_true&) // fallback scan angle +{ + using Point_3 = typename boost::property_traits::value_type; + using Vector_3 = typename boost::property_traits::value_type; + + Vector_3 direction; + Vector_3 vertical; + std::tie (direction, vertical) + = scanline_base (begin, end, point_map); + + // Estimate scanner position: + // average XY-projected point, located above + double mean_x = 0.; + double mean_y = 0.; + double max_z = -(std::numeric_limits::max)(); + std::size_t nb = 0; + + for (Iterator it = begin; it != end; ++ it) + { + const Point_3& p = get (point_map, *it); + mean_x += p.x(); + mean_y += p.y(); + max_z = (std::max)(max_z, p.z()); + ++ nb; + } + + Iterator last = end; -- last; + double length = CGAL::approximate_sqrt (CGAL::squared_distance + (get (point_map, *begin), + get (point_map, *last))); + + Point_3 scan_position (mean_x / nb, mean_y / nb, + max_z + length * 2); + + for (Iterator it = begin; it != end; ++ it) + { + Vector_3 line_of_sight (get(point_map, *it), scan_position); + const Vector_3& normal = get (normal_map, *it); + if (normal_along_scanline_is_inverted (normal, line_of_sight)) + put (normal_map, *it, -normal); + } +} + +} // namespace internal + +} // namespace Point_set_processing_3 +/// \endcond + +// ---------------------------------------------------------------------------- +// Public section +// ---------------------------------------------------------------------------- + +/** + \ingroup PkgPointSetProcessing3Algorithms + + orients the normals of the range of `points` by estimating a line + of sight and checking its consistency with the current normal orientation. + + \warning This function requires the input `points` to be ordered + along scanlines aligned on the XY-plane. It is typically designed + for 2.5D urban datasets acquired through, for example, airborne + LIDAR devices. + + First, scanlines are estimated as subranges of `points` by + iterating on `points`: + + - if the named parameter `scanline_id_map` is provided, the range + is cutted everytime the id changes. + + - if no scanline ID map is provided, a fallback method simply cuts + the range everytime 3 consecutive points form an acute angle on + the projected XY-plane. This fallback method gives suboptimal + results. + + Then, the line of sight (estimated vector between a point and the + position of the scanner at its time of acquisition) is estimated: + + - if `scan_angle` is provided, the line of sight can be directly + computed as a combination of the estimated scanline and of the + scan angle. + + - if no scan angle map is provided, then for each scanline, the + position of the scanner is estimated as being above of the + barycenter of the points of the scanline projected on the + XY-plane. This fallback method gives suboptimal results. + + Once the line of sight is estimated for each point, the normals are + oriented by checking, for each of them, if the line of sight and the + normal vector give a positive scalar product. If they don't, then + the normal vector is inverted. + + \note This method gives optimal results when `scanline_id_map` + and `scan_angle` are provided. Correct results may still be + produced in the absence of either one or both of these properties, + as long as the point set is ordered in 2.5D scanlines. + + \tparam PointRange is a model of `Range`. The value type of + its iterator is the key type of the named parameter `point_map`. + + \param points input point range. + \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below + + \cgalNamedParamsBegin + \cgalParamNBegin{point_map} + \cgalParamDescription{a property map associating points to the elements of the point set `points`} + \cgalParamType{a model of `ReadablePropertyMap` whose key type is the value type + of the iterator of `PointRange` and whose value type is `geom_traits::Point_3`} + \cgalParamDefault{`CGAL::Identity_property_map`} + \cgalParamNEnd + + \cgalParamNBegin{normal_map} + \cgalParamDescription{a property map associating normals to the + elements of the point set `points`} + \cgalParamType{a model of `WritablePropertyMap` whose key type + is the value type of the iterator of + `PointRange` and whose value type is + `geom_traits::Vector_3`} + \cgalParamNEnd + + \cgalParamNBegin{scan_angle_map} + \cgalParamDescription{a property map associating the angle of + acquisition (in degrees) to the elements of the point set + `points`} + \cgalParamType{a model of `ReadablePropertyMap` whose key type + is the value type of the iterator of + `PointRange` and whose value type is convertible + to `double`} + \cgalParamNEnd + + \cgalParamNBegin{scanline_id_map} + \cgalParamDescription{a property map associating a scanline ID + to the elements of the point set `points`. A scanline is + detected as a consecutive subrange of items in the input range + `point` whose ID are identical. IDs do not need to be unique, + they just need to be different for two consecutive + scanlines. The LAS property `scan_direction_flag` (whose values + are either 0 or 1 depending on the direction of the scanner) + can be used.} + \cgalParamType{a model of `ReadablePropertyMap` whose key type + is the value type of the iterator of + `PointRange` and whose value type is a model of + `EqualityComparable`} + \cgalParamNEnd + + \cgalParamNBegin{geom_traits} + \cgalParamDescription{an instance of a geometric traits class} + \cgalParamType{a model of `Kernel`} + \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`} + \cgalParamNEnd + \cgalNamedParamsEnd +*/ +template +void scanline_orient_normals (PointRange& points, const NamedParameters& np) +{ + using parameters::choose_parameter; + using parameters::get_parameter; + + using Iterator = typename PointRange::iterator; + using Value_type = typename std::iterator_traits::value_type; + using PointMap = typename CGAL::GetPointMap::type; + using NormalMap = typename Point_set_processing_3::GetNormalMap::type; + + using No_map = Constant_property_map; + + using ScanAngleMap = typename internal_np::Lookup_named_param_def + ::type; + using Fallback_scan_angle = Boolean_tag::value>; + + using ScanlineIDMap = typename internal_np::Lookup_named_param_def + ::type; + using Fallback_scanline_ID = Boolean_tag::value>; + + CGAL_static_assertion_msg(!(std::is_same::NoMap>::value), + "Error: no normal map"); + + PointMap point_map = choose_parameter(get_parameter(np, internal_np::point_map)); + NormalMap normal_map = choose_parameter(get_parameter(np, internal_np::normal_map)); + ScanAngleMap scan_angle_map = choose_parameter + (get_parameter(np, internal_np::scan_angle_map)); + ScanlineIDMap scanline_id_map = choose_parameter + (get_parameter(np, internal_np::scanline_id_map)); + + std::size_t nb_scanlines = 1; + +#ifdef CGAL_SCANORIENT_DUMP_RANDOM_SCANLINES + std::ofstream ofile ("scanlines.polylines.txt"); + ofile.precision(18); +#endif + + CGAL::Bbox_3 bbox = CGAL::bbox_3 + (boost::make_transform_iterator + (points.begin(), Property_map_to_unary_function(point_map)), + boost::make_transform_iterator + (points.end(), Property_map_to_unary_function(point_map))); + + double bbox_diagonal + = CGAL::approximate_sqrt((bbox.xmax() - bbox.xmin()) * (bbox.xmax() - bbox.xmin()) + + (bbox.ymax() - bbox.ymin()) * (bbox.ymax() - bbox.ymin()) + + (bbox.zmax() - bbox.zmin()) * (bbox.zmax() - bbox.zmin())); + double limit = 0.05 * bbox_diagonal; + + Iterator scanline_begin = points.begin(); + for (Iterator it = points.begin(); it != points.end(); ++ it) + { + bool force_cut = false; + if (it != points.begin()) + { + Iterator prev = it; -- prev; + force_cut = (CGAL::squared_distance (get (point_map, *prev), + get (point_map, *it)) > limit * limit); + } + + if (Point_set_processing_3::internal::is_end_of_scanline + (scanline_begin, it, point_map, scanline_id_map, + Fallback_scanline_ID()) || force_cut) + { + Point_set_processing_3::internal::orient_scanline + (scanline_begin, it, point_map, normal_map, + scan_angle_map, Fallback_scan_angle()); + +#ifdef CGAL_SCANORIENT_DUMP_RANDOM_SCANLINES + ofile << std::distance (scanline_begin, it); + for (Iterator it2 = scanline_begin; it2 != it; ++ it2) + ofile << " " << get (point_map, *it2); + ofile << std::endl; +#endif + + scanline_begin = it; + ++ nb_scanlines; + } + } + + Point_set_processing_3::internal::orient_scanline + (scanline_begin, points.end(), point_map, normal_map, + scan_angle_map, Fallback_scan_angle()); + +#ifdef CGAL_SCANLINE_ORIENT_VERBOSE + std::cerr << nb_scanlines << " scanline(s) identified (mean length = " + << std::size_t(points.size() / double(nb_scanlines)) + << " point(s))" << std::endl; +#endif +} + +template +void scanline_orient_normals (PointRange& points) +{ + return scanline_orient_normals (points, + CGAL::Point_set_processing_3::parameters::all_default(points)); +} + +} // namespace CGAL + +#endif // CGAL_EIGEN3_ENABLED + +#endif // CGAL_SCANLINE_ORIENT_NORMALS_H diff --git a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp index b4050ae9e92..f4168605e37 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -306,82 +307,156 @@ void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalOrientat if(points == NULL) return; - // Gets options - QMultipleInputDialog dialog ("Normal Orientation", mw); - QSpinBox* neighborhood = dialog.add ("Neighborhood Size: "); - neighborhood->setRange (1, 10000000); - neighborhood->setValue (18); + // Chose method + QMultipleInputDialog method ("Normal Orientation", mw); + QRadioButton* mst = method.add ("Orient by Minimum Spanning Tree"); + mst->setChecked(true); + QRadioButton* scanline = method.add ("Orient 2.5D airborne acquired scanlines"); + scanline->setChecked(false); - QRadioButton* use_seed_points = NULL; - QRadioButton* orient_selection = NULL; - - if (points->nb_selected_points() != 0) - { - use_seed_points = dialog.add ("Use selection as seed points and orient the unselected points"); - use_seed_points->setChecked(true); - orient_selection = dialog.add ("Orient selection"); - orient_selection->setChecked(false); - } - - if(!dialog.exec()) + if (!method.exec()) return; - QApplication::setOverrideCursor(Qt::BusyCursor); - QApplication::processEvents(); - - // First point to delete - Point_set::iterator first_unoriented_point = points->end(); - - //*************************************** - // normal orientation - //*************************************** - - CGAL::Timer task_timer; task_timer.start(); - std::cerr << "Orient normals with a Minimum Spanning Tree (k=" << neighborhood->value() << ")...\n"; - - // Tries to orient normals - if (points->nb_selected_points() != 0 && use_seed_points->isChecked()) + if (mst->isChecked()) { - std::vector constrained_map (points->size(), false); + // Gets options + QMultipleInputDialog dialog ("Normal Orientation", mw); + QSpinBox* neighborhood = dialog.add ("Neighborhood Size: "); + neighborhood->setRange (1, 10000000); + neighborhood->setValue (18); - for (Point_set::iterator it = points->first_selected(); it != points->end(); ++ it) - constrained_map[*it] = true; + QRadioButton* use_seed_points = NULL; + QRadioButton* orient_selection = NULL; - first_unoriented_point = - CGAL::mst_orient_normals(*points, - std::size_t(neighborhood->value()), - points->parameters(). - point_is_constrained_map(Vector_to_pmap(&constrained_map))); + if (points->nb_selected_points() != 0) + { + use_seed_points = dialog.add ("Use selection as seed points and orient the unselected points"); + use_seed_points->setChecked(true); + orient_selection = dialog.add ("Orient selection"); + orient_selection->setChecked(false); + } + + if(!dialog.exec()) + return; + + QApplication::setOverrideCursor(Qt::BusyCursor); + QApplication::processEvents(); + + // First point to delete + Point_set::iterator first_unoriented_point = points->end(); + + //*************************************** + // normal orientation + //*************************************** + + CGAL::Timer task_timer; task_timer.start(); + std::cerr << "Orient normals with a Minimum Spanning Tree (k=" << neighborhood->value() << ")...\n"; + + // Tries to orient normals + if (points->nb_selected_points() != 0 && use_seed_points->isChecked()) + { + std::vector constrained_map (points->size(), false); + + for (Point_set::iterator it = points->first_selected(); it != points->end(); ++ it) + constrained_map[*it] = true; + + first_unoriented_point = + CGAL::mst_orient_normals(*points, + std::size_t(neighborhood->value()), + points->parameters(). + point_is_constrained_map(Vector_to_pmap(&constrained_map))); + } + else + first_unoriented_point = + CGAL::mst_orient_normals(points->all_or_selection_if_not_empty(), + std::size_t(neighborhood->value()), + points->parameters()); + + std::size_t nb_unoriented_normals = std::distance(first_unoriented_point, points->end()); + std::size_t memory = CGAL::Memory_sizer().virtual_size(); + std::cerr << "Orient normals: " << nb_unoriented_normals << " point(s) with an unoriented normal are selected (" + << task_timer.time() << " seconds, " + << (memory>>20) << " Mb allocated)" + << std::endl; + + // Selects points with an unoriented normal + points->set_first_selected (first_unoriented_point); + + // Updates scene + item->invalidateOpenGLBuffers(); + scene->itemChanged(index); + + QApplication::restoreOverrideCursor(); + + // Warns user + if (nb_unoriented_normals > 0) + { + QMessageBox::information(NULL, + tr("Points with an unoriented normal"), + tr("%1 point(s) with an unoriented normal are selected.\nPlease orient them or remove them before running Poisson reconstruction.") + .arg(nb_unoriented_normals)); + } } - else - first_unoriented_point = - CGAL::mst_orient_normals(points->all_or_selection_if_not_empty(), - std::size_t(neighborhood->value()), - points->parameters()); - - std::size_t nb_unoriented_normals = std::distance(first_unoriented_point, points->end()); - std::size_t memory = CGAL::Memory_sizer().virtual_size(); - std::cerr << "Orient normals: " << nb_unoriented_normals << " point(s) with an unoriented normal are selected (" - << task_timer.time() << " seconds, " - << (memory>>20) << " Mb allocated)" - << std::endl; - - // Selects points with an unoriented normal - points->set_first_selected (first_unoriented_point); - - // Updates scene - item->invalidateOpenGLBuffers(); - scene->itemChanged(index); - - QApplication::restoreOverrideCursor(); - - // Warns user - if (nb_unoriented_normals > 0) + else // scanline method { - QMessageBox::information(NULL, - tr("Points with an unoriented normal"), - tr("%1 point(s) with an unoriented normal are selected.\nPlease orient them or remove them before running Poisson reconstruction.") - .arg(nb_unoriented_normals)); + QApplication::setOverrideCursor(Qt::BusyCursor); + QApplication::processEvents(); + + //*************************************** + // normal orientation + //*************************************** + + CGAL::Timer task_timer; task_timer.start(); + std::cerr << "Orient normals with along 2.5D scanlines..." << std::endl; + + Point_set::Property_map scan_angle; + Point_set::Property_map scan_direction_flag; + bool angle_found = false, flag_found = false; + + std::tie (scan_angle, angle_found) + = points->property_map("scan_angle"); + std::tie (scan_direction_flag, flag_found) + = points->property_map("scan_direction_flag"); + + if (!angle_found && !flag_found) + { + std::cerr << " using no additional properties" << std::endl; + CGAL::scanline_orient_normals(points->all_or_selection_if_not_empty(), + points->parameters()); + } + else if (!angle_found && flag_found) + { + std::cerr << " using scan direction flag" << std::endl; + CGAL::scanline_orient_normals(points->all_or_selection_if_not_empty(), + points->parameters(). + scanline_id_map (scan_direction_flag)); + } + else if (angle_found && !flag_found) + { + std::cerr << " using scan angle" << std::endl; + CGAL::scanline_orient_normals(points->all_or_selection_if_not_empty(), + points->parameters(). + scan_angle_map (scan_angle)); + } + else // if (angle_found && flag_found) + { + std::cerr << " using scan angle and direction flag" << std::endl; + CGAL::scanline_orient_normals(points->all_or_selection_if_not_empty(), + points->parameters(). + scan_angle_map (scan_angle). + scanline_id_map (scan_direction_flag)); + } + std::size_t memory = CGAL::Memory_sizer().virtual_size(); + std::cerr << "Orient normals: " + << task_timer.time() << " seconds, " + << (memory>>20) << " Mb allocated)" + << std::endl; + + // Updates scene + item->invalidateOpenGLBuffers(); + scene->itemChanged(index); + + QApplication::restoreOverrideCursor(); } } }