Merge pull request #5071 from sgiraudot/PSP-Scanline_orient_normals-GF

[Small Feature] Scanline Orient Normals

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)

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;

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
+following functions can be used:
-be used directly on input normals if their orientation is not
+
-consistent.
+- `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
 

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

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()`
 
 */
-

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 <I>Surface reconstruction from unorganized points</I> \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
 

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

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)
 

Point_set_processing_3/examples/Point_set_processing_3/orient_scanlines_example.cpp

@@ -0,0 +1,88 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/IO/read_las_points.h>
+#include <CGAL/IO/write_ply_points.h>
+#include <CGAL/jet_estimate_normals.h>
+#include <CGAL/scanline_orient_normals.h>
+
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point_3 = Kernel::Point_3;
+using Vector_3 = Kernel::Vector_3;
+using Point_with_info = std::tuple<Point_3, Vector_3, float, unsigned char>;
+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<Point_with_info>& 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<Point_with_info> 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<CGAL::Parallel_if_available_tag>
+    (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;
+}

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 <CGAL/license/Point_set_processing_3.h>
+
+#include <CGAL/squared_distance_3.h>
+
+#include <CGAL/boost/graph/Named_function_parameters.h>
+#include <CGAL/boost/graph/named_params_helper.h>
+
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+#include <boost/iterator/transform_iterator.hpp>
+
+#if defined(CGAL_EIGEN3_ENABLED) || defined(DOXYGEN_RUNNING)
+
+#include <Eigen/Dense>
+#include <Eigen/IterativeLinearSolvers>
+
+//#define CGAL_SCANLINE_ORIENT_VERBOSE
+
+namespace CGAL
+{
+
+/// \cond SKIP_IN_MANUAL
+namespace Point_set_processing_3
+{
+
+namespace internal
+{
+
+template <typename Vector_3>
+const Vector_3& vertical_vector()
+{
+  static Vector_3 v(0, 0, 1);
+  return v;
+}
+
+template <typename Iterator, typename PointMap,
+          typename ScanlineIDMap>
+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 <typename Iterator, typename PointMap,
+          typename ScanlineIDMap>
+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<PointMap>::value_type;
+  using Vector_3 = typename Kernel_traits<Point_3>::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 <typename Iterator, typename PointMap>
+std::pair<typename Kernel_traits<typename boost::property_traits
+                                 <PointMap>::value_type>::Kernel::Vector_3,
+          typename Kernel_traits<typename boost::property_traits
+                                 <PointMap>::value_type>::Kernel::Vector_3>
+scanline_base (Iterator begin, Iterator end,
+               PointMap point_map)
+{
+  using Point_3 = typename boost::property_traits<PointMap>::value_type;
+  using Kernel = typename Kernel_traits<Point_3>::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<PointMap>(point_map)),
+      boost::make_transform_iterator
+     (end, Property_map_to_unary_function<PointMap>(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<PointMap>(point_map)),
+     boost::make_transform_iterator
+     (end, Property_map_to_unary_function<PointMap>(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<Vector_3>() < 0)
+    vertical = -vertical;
+
+  vertical = vertical / CGAL::approximate_sqrt(vertical.squared_length());
+
+  if (std::acos(vertical * vertical_vector<Vector_3>()) < 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<Vector_3>());
+}
+
+template <typename Vector_3>
+bool normal_along_scanline_is_inverted
+(const Vector_3& normal, const Vector_3& line_of_sight)
+{
+  return (line_of_sight * normal < 0);
+}
+
+template <typename Iterator, typename PointMap, typename Vector_3>
+std::pair<typename boost::property_traits<PointMap>::value_type, bool>
+estimate_scan_position (Iterator begin, Iterator end, PointMap point_map,
+                        const std::vector<Vector_3>& lines_of_sight)
+{
+  using Point_3 = typename boost::property_traits<PointMap>::value_type;
+  typedef Eigen::Matrix3d Matrix;
+  typedef Eigen::Vector3d Vector;
+  typedef Eigen::ConjugateGradient<Matrix> 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 <typename Iterator, typename PointMap, typename NormalMap,
+          typename ScanAngleMap>
+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<PointMap>::value_type;
+  using Vector_3 = typename boost::property_traits<NormalMap>::value_type;
+
+  Vector_3 direction;
+  Vector_3 vertical;
+  std::tie (direction, vertical)
+    = scanline_base (begin, end, point_map);
+
+  std::vector<Vector_3> 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<double>(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<double>(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 <typename Iterator, typename PointMap, typename NormalMap,
+          typename ScanAngleMap>
+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<PointMap>::value_type;
+  using Vector_3 = typename boost::property_traits<NormalMap>::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<double>::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<geom_traits::Point_3>`}
+     \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 <typename PointRange, typename NamedParameters>
+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<Iterator>::value_type;
+  using PointMap = typename CGAL::GetPointMap<PointRange, NamedParameters>::type;
+  using NormalMap = typename Point_set_processing_3::GetNormalMap<PointRange, NamedParameters>::type;
+
+  using No_map = Constant_property_map<Value_type, int>;
+
+  using ScanAngleMap = typename internal_np::Lookup_named_param_def
+    <internal_np::scan_angle_t, NamedParameters, No_map>::type;
+  using Fallback_scan_angle = Boolean_tag<std::is_same<ScanAngleMap, No_map>::value>;
+
+  using ScanlineIDMap = typename internal_np::Lookup_named_param_def
+    <internal_np::scanline_id_t, NamedParameters, No_map>::type;
+  using Fallback_scanline_ID = Boolean_tag<std::is_same<ScanlineIDMap, No_map>::value>;
+
+  CGAL_static_assertion_msg(!(std::is_same<NormalMap,
+                              typename Point_set_processing_3::GetNormalMap
+                              <PointRange, NamedParameters>::NoMap>::value),
+                            "Error: no normal map");
+
+  PointMap point_map = choose_parameter<PointMap>(get_parameter(np, internal_np::point_map));
+  NormalMap normal_map = choose_parameter<NormalMap>(get_parameter(np, internal_np::normal_map));
+  ScanAngleMap scan_angle_map = choose_parameter<ScanAngleMap>
+    (get_parameter(np, internal_np::scan_angle_map));
+  ScanlineIDMap scanline_id_map = choose_parameter<ScanlineIDMap>
+    (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<PointMap>(point_map)),
+     boost::make_transform_iterator
+     (points.end(), Property_map_to_unary_function<PointMap>(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 <typename PointRange>
+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

Polyhedron/demo/Polyhedron/Plugins/Point_set/Point_set_normal_estimation_plugin.cpp

@@ -7,6 +7,7 @@
 #include <CGAL/jet_estimate_normals.h>
 #include <CGAL/vcm_estimate_normals.h>
 #include <CGAL/mst_orient_normals.h>
+#include <CGAL/scanline_orient_normals.h>
 #include <CGAL/Timer.h>
 #include <CGAL/Memory_sizer.h>
 
@@ -306,6 +307,18 @@ void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalOrientat
     if(points == NULL)
         return;
 
+    // Chose method
+    QMultipleInputDialog method ("Normal Orientation", mw);
+    QRadioButton* mst = method.add<QRadioButton> ("Orient by Minimum Spanning Tree");
+    mst->setChecked(true);
+    QRadioButton* scanline = method.add<QRadioButton> ("Orient 2.5D airborne acquired scanlines");
+    scanline->setChecked(false);
+
+    if (!method.exec())
+      return;
+
+    if (mst->isChecked())
+    {
       // Gets options
       QMultipleInputDialog dialog ("Normal Orientation", mw);
       QSpinBox* neighborhood = dialog.add<QSpinBox> ("Neighborhood Size: ");
@@ -384,6 +397,68 @@ void Polyhedron_demo_point_set_normal_estimation_plugin::on_actionNormalOrientat
                                  .arg(nb_unoriented_normals));
       }
     }
+    else // scanline method
+    {
+      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<float> scan_angle;
+      Point_set::Property_map<unsigned char> scan_direction_flag;
+      bool angle_found = false, flag_found = false;
+
+      std::tie (scan_angle, angle_found)
+        = points->property_map<float>("scan_angle");
+      std::tie (scan_direction_flag, flag_found)
+        = points->property_map<unsigned char>("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();
+    }
+  }
 }