1 // Copyright (C) 2013 INRIA - Sophia Antipolis (France). 2 // Copyright (c) 2017 GeometryFactory Sarl (France). 3 // 4 // This file is part of CGAL (www.cgal.org). 5 // 6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h $ 7 // $Id: Scale_space_surface_reconstruction_3.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial 9 // 10 // Author(s): Thijs van Lankveld, Simon Giraudot 11 12 13 #ifndef CGAL_SCALE_SPACE_SURFACE_RECONSTRUCTION_3_H 14 #define CGAL_SCALE_SPACE_SURFACE_RECONSTRUCTION_3_H 15 16 #include <CGAL/license/Scale_space_reconstruction_3.h> 17 18 #include <CGAL/Scale_space_reconstruction_3/Weighted_PCA_smoother.h> 19 #include <CGAL/Scale_space_reconstruction_3/Alpha_shape_mesher.h> 20 21 22 namespace CGAL 23 { 24 25 26 /// computes a triangulated surface mesh interpolating a point set. 27 /** \ingroup PkgScaleSpaceReconstruction3Classes 28 * 29 * Scale space reconstruction consists in combining a smoothing 30 * algorithm (see `Scale_space_reconstruction_3::Smoother`) that is 31 * used to generate a user-specified number of scales with a meshing 32 * algorithm (see `Scale_space_reconstruction_3::Mesher`). 33 * 34 * The class stores the point set at the current scale and the 35 * reconstructed surface. User can either use the reconstructed 36 * surface at the current scale (smoothed version) or propagate the 37 * connectivity of the facets back to the original point set. This 38 * latest option makes it possible to reconstruct an interpolative 39 * surface on a noisy point set. 40 * 41 * \tparam Geom_traits is the geometric traits class. It must be a 42 * model of `DelaunayTriangulationTraits_3`. It must have a 43 * `RealEmbeddable` field number type. Generally, 44 * `Exact_predicates_inexact_constructions_kernel` is preferred. 45 */ 46 template <typename Geom_traits> 47 class Scale_space_surface_reconstruction_3 48 { 49 public: 50 51 typedef typename Geom_traits::FT FT; ///< defines the field number type. 52 typedef typename Geom_traits::Point_3 Point; ///< defines the point type. 53 typedef std::array<std::size_t, 3> Facet; ///< defines a facet of the surface (triple of point indices). 54 55 #ifdef DOXYGEN_RUNNING 56 typedef unspecified_type Point_iterator; ///< defines an iterator over the points. 57 typedef const unspecified_type Point_const_iterator; ///< defines a constant iterator over the points. 58 #else 59 typedef typename std::vector<Point> Point_vector; 60 typedef typename Point_vector::iterator Point_iterator; 61 typedef typename Point_vector::const_iterator Point_const_iterator; 62 #endif 63 64 #ifdef DOXYGEN_RUNNING 65 typedef unspecified_type Facet_iterator; ///< defines an iterator over the points. 66 typedef const unspecified_type Facet_const_iterator; ///< defines a constant iterator over the points. 67 #else 68 typedef typename std::vector<Facet> Facet_vector; 69 typedef typename Facet_vector::iterator Facet_iterator; 70 typedef typename Facet_vector::const_iterator Facet_const_iterator; 71 #endif 72 73 // Default algorithms used (same as in old API) 74 typedef Scale_space_reconstruction_3::Weighted_PCA_smoother<Geom_traits> Weighted_PCA_smoother; 75 typedef Scale_space_reconstruction_3::Alpha_shape_mesher<Geom_traits> Alpha_shape_mesher; 76 77 private: 78 79 Point_vector m_points; 80 Facet_vector m_facets; 81 82 FT m_internal_squared_radius; // For backward compatibility 83 84 public: 85 86 /// @{ 87 /// \name Initialization 88 89 /** 90 * Empty constructor. 91 */ Scale_space_surface_reconstruction_3()92 Scale_space_surface_reconstruction_3 () : m_internal_squared_radius(0.) { } 93 94 /** 95 * Constructs a reconstruction object with the given point range. 96 * 97 * \note This constructor is equivalent to constructing an empty 98 * object and then calling <code>[insert(begin, end)](\ref insert)</code>. 99 * 100 * \tparam InputIterator is an iterator over the point collection. 101 * The value type of the iterator must be a `Point`. 102 * 103 * \param begin is an iterator to the first point of the collection. 104 * \param end is a past-the-end iterator for the point collection. 105 */ 106 template <typename InputIterator> Scale_space_surface_reconstruction_3(InputIterator begin,InputIterator end)107 Scale_space_surface_reconstruction_3 (InputIterator begin, InputIterator end) 108 : m_internal_squared_radius (0.) 109 { 110 insert (begin, end); 111 } 112 113 /// inserts a point into the scale-space at the current scale. 114 /** \param p is the point to insert. 115 * 116 * \note Inserting the point does not automatically construct or 117 * update the surface. 118 * 119 * \note In order to construct the surface, call 120 * `#reconstruct_surface()`. 121 * 122 * \sa `insert(InputIterator begin, InputIterator end)`. 123 */ insert(const Point & p)124 void insert (const Point& p) 125 { 126 m_points.push_back (p); 127 } 128 129 /// inserts a collection of points into the scale-space at the current scale. 130 /** \tparam InputIterator is an iterator over the point collection. 131 * The value type of the iterator must be a `Point`. 132 * 133 * \param begin is an iterator to the first point of the collection. 134 * \param end is a past-the-end iterator for the point collection. 135 * 136 * \note Inserting the points does not automatically construct or 137 * update the surface. 138 * 139 * \note In order to construct the surface, call 140 * `reconstruct_surface()` after inserting the points. 141 * 142 * \sa `insert(const Point& p)`. 143 */ 144 template <typename InputIterator> insert(InputIterator begin,InputIterator end)145 void insert (InputIterator begin, InputIterator end) 146 { 147 std::copy (begin, end, std::back_inserter (m_points)); 148 } 149 150 /// @} 151 152 /// @{ 153 /// \name Algorithms 154 155 /// increases the scale by a number of iterations. 156 /** Each iteration the scale is increased, the points set at a higher scale 157 * is computed. At a higher scale, the points set is smoother. 158 * 159 * \note If no smoothing algorithm is specified, the default 160 * `Weighted_PCA_smoother` is used. 161 * 162 * \tparam Smoother model of `Scale_space_reconstruction_3::Smoother`. 163 * 164 * \param iterations is the number of iterations to perform. If 165 * `iterations` is 0, nothing happens. 166 * \param smoother the smoother object. 167 * 168 * \note This method processes the point set at the current scale. The 169 * points can be set with <code>[insert(begin, end)](\ref insert)</code>. 170 * 171 * \note If the surface was already constructed, increasing the scale 172 * will not automatically adjust the surface. 173 * 174 * \sa `reconstruct_surface()`. 175 */ 176 template <typename Smoother> 177 void increase_scale (std::size_t iterations = 1, const Smoother& smoother = Smoother()) 178 { 179 for (std::size_t i = 0; i < iterations; ++ i) 180 const_cast<Smoother&>(smoother) (m_points.begin(), m_points.end()); 181 } 182 183 /// \cond SKIP_IN_MANUAL 184 void increase_scale (std::size_t iterations = 1) 185 { 186 Weighted_PCA_smoother smoother; 187 increase_scale (iterations, smoother); 188 m_internal_squared_radius = smoother.squared_radius(); 189 } 190 /// \endcond 191 192 /// constructs a triangle mesh from the point set at a fixed scale. 193 /** The order of the points at the current scale is the same as the order 194 * at the original scale, meaning that the surface can interpolate the 195 * point set at the original scale by applying the indices of the surface 196 * triangles to the original point set. 197 * 198 * After construction, the facets of the surface can be iterated using 199 * `facets_begin()` and `facets_end()`. 200 * 201 * \note If no meshing algorithm is specified, the default 202 * `Alpha_shape_mesher` is used. 203 * 204 * \tparam Mesher model of `Scale_space_reconstruction_3::Mesher`. 205 * 206 * \param mesher the mesher object. 207 * 208 * \note This method processes the point set at the current scale. The 209 * points can be set with <code>[insert(begin, end)](\ref insert)</code>. 210 * 211 * \sa `increase_scale()`. 212 */ 213 template <typename Mesher> 214 void reconstruct_surface (const Mesher& mesher = Mesher()) 215 { 216 m_facets.clear(); 217 const_cast<Mesher&>(mesher) (m_points.begin(), m_points.end(), std::back_inserter (m_facets)); 218 } 219 220 /// \cond SKIP_IN_MANUAL reconstruct_surface()221 void reconstruct_surface () 222 { 223 CGAL_assertion (m_internal_squared_radius != 0.); 224 reconstruct_surface (Alpha_shape_mesher(m_internal_squared_radius)); 225 } 226 /// \endcond 227 228 /// @} 229 230 /// @{ 231 /// \name Output 232 233 /// gives the number of points of the surface. number_of_points()234 std::size_t number_of_points() const { return m_points.size(); } 235 236 /// gives an iterator to the first point at the current scale. 237 /** \warning Changes to the scale-space do not cause an automatic update to 238 * the surface. 239 */ points_begin()240 Point_iterator points_begin() { return m_points.begin(); } 241 242 /// gives a past-the-end iterator of the points at the current scale. 243 /** \warning Changes to the scale-space do not cause an automatic update to 244 * the surface. 245 */ points_end()246 Point_iterator points_end() { return m_points.end(); } 247 248 /// gives an iterator to the first point at the current scale. points_begin()249 Point_const_iterator points_begin() const { return m_points.begin(); } 250 251 /// gives a past-the-end iterator of the points at the current scale. points_end()252 Point_const_iterator points_end() const { return m_points.end(); } 253 254 /// gives the number of facets of the surface. number_of_facets()255 std::size_t number_of_facets() const { return m_facets.size(); } 256 257 /// gives an iterator to the first triple in the surface. 258 /** \warning Changes to the surface may change its topology. 259 */ facets_begin()260 Facet_iterator facets_begin() { return m_facets.begin(); } 261 262 /// gives a past-the-end iterator of the triples in the surface. 263 /** \warning Changes to the surface may change its topology. 264 */ facets_end()265 Facet_iterator facets_end() { return m_facets.end(); } 266 267 /// gives an iterator to the first triple in the surface. facets_begin()268 Facet_const_iterator facets_begin() const { return m_facets.begin(); } 269 270 /// gives a past-the-end iterator of the triples in the surface. facets_end()271 Facet_const_iterator facets_end() const { return m_facets.end(); } 272 273 /// @} 274 }; 275 276 277 template <typename Geom_traits> 278 std::ostream& operator<< (std::ostream& os, const CGAL::Scale_space_surface_reconstruction_3<Geom_traits>& ss) 279 { 280 typedef CGAL::Scale_space_surface_reconstruction_3<Geom_traits> Reconstruction; 281 282 os << "OFF" << std::endl 283 << ss.number_of_points() << " " << ss.number_of_facets() << " 0" << std::endl; 284 for (typename Reconstruction::Point_const_iterator it = ss.points_begin(); 285 it != ss.points_end(); ++ it) 286 os << *it << std::endl; 287 for (typename Reconstruction::Facet_const_iterator it = ss.facets_begin(); 288 it != ss.facets_end(); ++ it) 289 os << "3 " << (*it)[0] << " " << (*it)[1] << " " << (*it)[2] << std::endl; 290 return os; 291 } 292 293 } // namespace CGAL 294 295 template< typename T > 296 std::ostream& 297 operator<<( std::ostream& os, const std::array< T, 3 >& t ) { 298 return os << t[0] << " " << t[1] << " " << t[2]; 299 } 300 301 template< typename T > 302 std::istream& 303 operator>>( std::istream& is, std::array< T, 3 >& t ) { 304 return is >> get<0>(t) >> get<1>(t) >> get<2>(t); 305 } 306 307 #endif // CGAL_SCALE_SPACE_SURFACE_RECONSTRUCTION_3_H 308