1 // Copyright (c) 2018 Liangliang Nan. All rights reserved. 2 // 3 // This file is part of CGAL (www.cgal.org) 4 // 5 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Polygonal_surface_reconstruction/include/CGAL/Polygonal_surface_reconstruction.h $ 6 // $Id: Polygonal_surface_reconstruction.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot 7 // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial 8 // 9 // Author(s) : Liangliang Nan 10 11 #ifndef CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_H 12 #define CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_H 13 14 #include <CGAL/license/Polygonal_surface_reconstruction.h> 15 16 #include <CGAL/bounding_box.h> 17 #include <CGAL/property_map.h> 18 #include <CGAL/internal/hypothesis.h> 19 #include <CGAL/internal/compute_confidences.h> 20 #include <CGAL/internal/point_set_with_planes.h> 21 22 #include <unordered_map> 23 24 /*! 25 \file Polygonal_surface_reconstruction.h 26 */ 27 28 namespace CGAL { 29 30 /*! 31 \ingroup PkgPolygonalSurfaceReconstruction 32 33 \brief 34 35 Implementation of the Polygonal Surface Reconstruction method. 36 37 Given a set of 3D points with unoriented normals sampled 38 from the outer boundary of a piecewise planar object, the Polygonal Surface 39 Reconstruction method \cgalCite{nan2017polyfit} outputs a simplified and 40 watertight surface mesh interpolating the input point set. 41 42 The method first generates a set of face candidates by intersecting the planar 43 primitives. Then an optimal subset of the candidate faces is selected through 44 optimization under hard constraints that enforce the final model to be manifold 45 and watertight. 46 47 The reconstruction assumes the planar segmentation of the point cloud is 48 provided in the input. 49 50 \tparam GeomTraits a geometric traits class, model of Kernel 51 */ 52 template <class GeomTraits> 53 class Polygonal_surface_reconstruction 54 { 55 public: 56 57 /// \name Types 58 59 typedef typename GeomTraits::FT FT; ///< number type. 60 typedef typename GeomTraits::Point_3 Point; ///< point type. 61 typedef typename GeomTraits::Vector_3 Vector; ///< vector type. 62 typedef typename GeomTraits::Plane_3 Plane; ///< plane type. 63 64 private: 65 66 // Polygon mesh storing the candidate faces 67 typedef CGAL::Surface_mesh<Point> Polygon_mesh; 68 69 typedef typename Polygon_mesh::Face_index Face_descriptor; 70 typedef typename Polygon_mesh::Vertex_index Vertex_descriptor; 71 typedef typename Polygon_mesh::Edge_index Edge_descriptor; 72 typedef typename Polygon_mesh::Halfedge_index Halfedge_descriptor; 73 74 // Public methods 75 public: 76 77 /// \name Creation 78 79 /*! 80 Creates a Polygonal Surface Reconstruction object. 81 After construction, candidate faces are generated and point/face confidence values are 82 computed, allowing to reuse them in the subsequent reconstruction step with different parameters. 83 84 \tparam PointRange is the range of input points, model of `ConstRange`. 85 \tparam PointMap is a model of `ReadablePropertyMap` with value type `GeomTraits::Point_3`. 86 \tparam NormalMap is a model of `ReadablePropertyMap` with value type `GeomTraits::Vector_3`. 87 \tparam IndexMap is a model of `ReadablePropertyMap` with value type `int`. 88 89 \param points range of input points. 90 \param point_map property map: value_type of `typename PointRange::const_iterator` -> `Point_3` 91 \param normal_map property map: value_type of `typename PointRange::const_iterator` -> `Vector_3` 92 \param index_map property map: value_type of `typename PointRange::const_iterator` -> `int`, 93 denoting the index of the plane it belongs to (-1 if the point is not assigned to a plane) 94 */ 95 template < 96 typename PointRange, 97 typename PointMap, 98 typename NormalMap, 99 typename IndexMap 100 > 101 Polygonal_surface_reconstruction( 102 const PointRange& points, 103 PointMap point_map, 104 NormalMap normal_map, 105 IndexMap index_map 106 ); 107 108 109 /// \name Operations 110 111 /** Reconstructs a watertight polygonal mesh model. 112 113 \tparam MixedIntegerProgramTraits a model of `MixedIntegerProgramTraits` 114 115 \tparam PolygonMesh a model of `MutableFaceGraph` 116 117 \return `true` if the reconstruction succeeded, `false` otherwise. 118 */ 119 template <typename MixedIntegerProgramTraits, typename PolygonMesh> 120 bool reconstruct( 121 PolygonMesh& output_mesh, ///< the final reconstruction result 122 double wt_fitting = 0.43, ///< weight for the data fitting term. 123 double wt_coverage = 0.27, ///< weight for the point coverage term. 124 double wt_complexity = 0.30 ///< weight for the model complexity term. 125 ); 126 127 /*! Gives the user the possibility to access the intermediate candidate faces 128 (i.e., the faces induced by the intersection of the supporting planes). 129 \tparam PolygonMesh a model of `MutableFaceGraph`. 130 */ 131 template <typename PolygonMesh> 132 void output_candidate_faces(PolygonMesh& candidate_faces) const; 133 134 /// Gets the error message (if reconstruction failed). error_message()135 const std::string& error_message() const { return error_message_; } 136 137 // Data members. 138 private: 139 internal::Hypothesis<GeomTraits> hypothesis_; 140 141 // The generated candidate faces stored as a polygon mesh 142 Polygon_mesh candidate_faces_; 143 144 std::string error_message_; 145 146 private: // Copying is not allowed 147 Polygonal_surface_reconstruction(const Polygonal_surface_reconstruction& psr); 148 149 }; // end of Polygonal_surface_reconstruction 150 151 152 //////////////////////////////////////////////////////////////////////////s 153 154 // implementations 155 156 template <class GeomTraits> 157 158 template < 159 typename PointRange, 160 typename PointMap, 161 typename NormalMap, 162 typename IndexMap 163 > Polygonal_surface_reconstruction(const PointRange & points,PointMap point_map,NormalMap normal_map,IndexMap index_map)164 Polygonal_surface_reconstruction<GeomTraits>::Polygonal_surface_reconstruction( 165 const PointRange& points, 166 PointMap point_map, 167 NormalMap normal_map, 168 IndexMap index_map 169 ) : error_message_("") 170 { 171 if (points.empty()) { 172 error_message_ = "empty input points"; 173 return; 174 } 175 176 typedef internal::Planar_segment<GeomTraits> Planar_segment; 177 typedef internal::Point_set_with_planes<GeomTraits> Point_set_with_planes; 178 179 Point_set_with_planes point_set(points, point_map, normal_map, index_map); 180 181 const std::vector< Planar_segment* >& planar_segments = point_set.planar_segments(); 182 if (planar_segments.size() < 4) { 183 error_message_ = "at least 4 planes required to reconstruct a closed surface mesh (only " 184 + std::to_string(planar_segments.size()) + " provided)"; 185 return; 186 } 187 188 hypothesis_.generate(point_set, candidate_faces_); 189 190 typedef internal::Candidate_confidences<GeomTraits> Candidate_confidences; 191 Candidate_confidences conf; 192 conf.compute(point_set, candidate_faces_); 193 } 194 195 196 template <class GeomTraits> 197 template <typename PolygonMesh> output_candidate_faces(PolygonMesh & candidate_faces)198 void Polygonal_surface_reconstruction<GeomTraits>::output_candidate_faces(PolygonMesh& candidate_faces) const { 199 candidate_faces.clear(); // make sure it is empty. 200 CGAL::copy_face_graph(candidate_faces_, candidate_faces); 201 } 202 203 204 template <class GeomTraits> 205 template <typename MixedIntegerProgramTraits, typename PolygonMesh> reconstruct(PolygonMesh & output_mesh,double wt_fitting,double wt_coverage,double wt_complexity)206 bool Polygonal_surface_reconstruction<GeomTraits>::reconstruct( 207 PolygonMesh& output_mesh, 208 double wt_fitting /* = 0.43 */, 209 double wt_coverage /* = 0.27 */, 210 double wt_complexity /* = 0.30 */) 211 { 212 if (!error_message_.empty()) { // an error has occurred in the constructor 213 return false; 214 } 215 216 if (candidate_faces_.num_faces() < 4) { 217 error_message_ = "at least 4 candidate faces required to reconstruct a closed surface mesh (only " 218 + std::to_string(candidate_faces_.num_faces()) + " computed)"; 219 return false; 220 } 221 222 typedef typename internal::Hypothesis<GeomTraits>::Adjacency Adjacency; 223 const Adjacency& adjacency = hypothesis_.extract_adjacency(candidate_faces_); 224 225 // Internal data structure 226 Polygon_mesh target_mesh = candidate_faces_; 227 228 // The number of supporting points of each face 229 typename Polygon_mesh::template Property_map<Face_descriptor, std::size_t> face_num_supporting_points = 230 target_mesh.template add_property_map<Face_descriptor, std::size_t>("f:num_supporting_points").first; 231 232 // The area of each face 233 typename Polygon_mesh::template Property_map<Face_descriptor, FT> face_areas = 234 target_mesh.template add_property_map<Face_descriptor, FT>("f:face_area").first; 235 236 // The point covered area of each face 237 typename Polygon_mesh::template Property_map<Face_descriptor, FT> face_covered_areas = 238 target_mesh.template add_property_map<Face_descriptor, FT>("f:covered_area").first; 239 240 // The supporting plane of each face 241 typename Polygon_mesh::template Property_map<Face_descriptor, const Plane*> face_supporting_planes = 242 target_mesh.template add_property_map<Face_descriptor, const Plane*>("f:supp_plane").first; 243 244 // Gives each face an index 245 typename Polygon_mesh::template Property_map<Face_descriptor, std::size_t> face_indices = 246 target_mesh.template add_property_map<Face_descriptor, std::size_t>("f:index").first; 247 248 double total_points = 0.0; 249 std::size_t idx = 0; 250 for(auto f : target_mesh.faces()) { 251 total_points += face_num_supporting_points[f]; 252 face_indices[f] = idx; 253 ++idx; 254 } 255 256 257 typedef MixedIntegerProgramTraits MIP_Solver; 258 typedef typename MixedIntegerProgramTraits::Variable Variable; 259 typedef typename MixedIntegerProgramTraits::Linear_objective Linear_objective; 260 typedef typename MixedIntegerProgramTraits::Linear_constraint Linear_constraint; 261 262 MIP_Solver solver; 263 264 // Adds variables 265 266 // Binary variables: 267 // x[0] ... x[num_faces - 1] : binary labels of all the input faces 268 // x[num_faces] ... x[num_faces + num_edges - 1] : binary labels of all the intersecting edges (remain or not) 269 // x[num_faces + num_edges] ... x[num_faces + num_edges + num_edges] : binary labels of corner edges (sharp edge of not) 270 271 std::size_t num_faces = target_mesh.number_of_faces(); 272 std::size_t num_edges(0); 273 274 typedef typename internal::Hypothesis<GeomTraits>::Intersection Intersection; 275 276 std::unordered_map<const Intersection*, std::size_t> edge_usage_status; // keep or remove an intersecting edges 277 for (std::size_t i = 0; i < adjacency.size(); ++i) { 278 const Intersection& fan = adjacency[i]; 279 if (fan.size() == 4) { 280 std::size_t var_idx = num_faces + num_edges; 281 edge_usage_status[&fan] = var_idx; 282 ++num_edges; 283 } 284 } 285 286 std::size_t total_variables = num_faces + num_edges + num_edges; 287 288 const std::vector<Variable*>& variables = solver.create_variables(total_variables); 289 for (std::size_t i = 0; i < total_variables; ++i) { 290 Variable* v = variables[i]; 291 v->set_variable_type(Variable::BINARY); 292 } 293 294 // Adds objective 295 296 const typename Polygon_mesh::template Property_map<Vertex_descriptor, Point>& coords = target_mesh.points(); 297 std::vector<Point> vertices(target_mesh.number_of_vertices()); 298 idx = 0; 299 for(auto v : target_mesh.vertices()) { 300 vertices[idx] = coords[v]; 301 ++idx; 302 } 303 304 typedef typename GeomTraits::Iso_cuboid_3 Box; 305 306 const Box& box = CGAL::bounding_box(vertices.begin(), vertices.end()); 307 FT dx = box.xmax() - box.xmin(); 308 FT dy = box.ymax() - box.ymin(); 309 FT dz = box.zmax() - box.zmin(); 310 FT box_area = FT(2.0) * (dx * dy + dy * dz + dz * dx); 311 312 // Chooses a better scale: all actual values multiplied by total number of points 313 double coeff_data_fitting = wt_fitting; 314 double coeff_coverage = total_points * wt_coverage / box_area; 315 double coeff_complexity = total_points * wt_complexity / double(adjacency.size()); 316 317 Linear_objective * objective = solver.create_objective(Linear_objective::MINIMIZE); 318 319 std::unordered_map<const Intersection*, std::size_t> edge_sharp_status; // the edge is sharp or not 320 std::size_t num_sharp_edges = 0; 321 for (std::size_t i = 0; i < adjacency.size(); ++i) { 322 const Intersection& fan = adjacency[i]; 323 if (fan.size() == 4) { 324 std::size_t var_idx = num_faces + num_edges + num_sharp_edges; 325 edge_sharp_status[&fan] = var_idx; 326 327 // Accumulates model complexity term 328 objective->add_coefficient(variables[var_idx], coeff_complexity); 329 ++num_sharp_edges; 330 } 331 } 332 CGAL_assertion(num_edges == num_sharp_edges); 333 334 for(auto f : target_mesh.faces()) { 335 std::size_t var_idx = face_indices[f]; 336 337 // Accumulates data fitting term 338 std::size_t num = face_num_supporting_points[f]; 339 objective->add_coefficient(variables[var_idx], -coeff_data_fitting * num); 340 341 // Accumulates model coverage term 342 double uncovered_area = (face_areas[f] - face_covered_areas[f]); 343 objective->add_coefficient(variables[var_idx], coeff_coverage * uncovered_area); 344 } 345 346 // Adds constraints: the number of faces associated with an edge must be either 2 or 0 347 std::size_t var_edge_used_idx = 0; 348 for (std::size_t i = 0; i < adjacency.size(); ++i) { 349 Linear_constraint* c = solver.create_constraint(0.0, 0.0); 350 const Intersection& fan = adjacency[i]; 351 for (std::size_t j = 0; j < fan.size(); ++j) { 352 Face_descriptor f = target_mesh.face(fan[j]); 353 std::size_t var_idx = face_indices[f]; 354 c->add_coefficient(variables[var_idx], 1.0); 355 } 356 357 if (fan.size() == 4) { 358 std::size_t var_idx = num_faces + var_edge_used_idx; 359 c->add_coefficient(variables[var_idx], -2.0); // 360 ++var_edge_used_idx; 361 } 362 else { // boundary edge 363 // will be set to 0 (i.e., we don't allow open surface) 364 } 365 } 366 367 // Adds constraints: for the sharp edges. The explanation of posing this constraint can be found here: 368 // https://user-images.githubusercontent.com/15526536/30185644-12085a9c-942b-11e7-831d-290dd2a4d50c.png 369 double M = 1.0; 370 for (std::size_t i = 0; i < adjacency.size(); ++i) { 371 const Intersection& fan = adjacency[i]; 372 if (fan.size() != 4) 373 continue; 374 375 // If an edge is sharp, the edge must be selected first: 376 // X[var_edge_usage_idx] >= X[var_edge_sharp_idx] 377 Linear_constraint* c = solver.create_constraint(0.0); 378 std::size_t var_edge_usage_idx = edge_usage_status[&fan]; 379 c->add_coefficient(variables[var_edge_usage_idx], 1.0); 380 std::size_t var_edge_sharp_idx = edge_sharp_status[&fan]; 381 c->add_coefficient(variables[var_edge_sharp_idx], -1.0); 382 383 for (std::size_t j = 0; j < fan.size(); ++j) { 384 Face_descriptor f1 = target_mesh.face(fan[j]); 385 const Plane* plane1 = face_supporting_planes[f1]; 386 std::size_t fid1 = face_indices[f1]; 387 for (std::size_t k = j + 1; k < fan.size(); ++k) { 388 Face_descriptor f2 = target_mesh.face(fan[k]); 389 const Plane* plane2 = face_supporting_planes[f2]; 390 std::size_t fid2 = face_indices[f2]; 391 392 if (plane1 != plane2) { 393 // The constraint is: 394 //X[var_edge_sharp_idx] + M * (3 - (X[fid1] + X[fid2] + X[var_edge_usage_idx])) >= 1 395 // which equals to 396 //X[var_edge_sharp_idx] - M * X[fid1] - M * X[fid2] - M * X[var_edge_usage_idx] >= 1 - 3M 397 c = solver.create_constraint(1.0 - 3.0 * M); 398 c->add_coefficient(variables[var_edge_sharp_idx], 1.0); 399 c->add_coefficient(variables[fid1], -M); 400 c->add_coefficient(variables[fid2], -M); 401 c->add_coefficient(variables[var_edge_usage_idx], -M); 402 } 403 } 404 } 405 } 406 407 // Optimization 408 409 if (solver.solve()) { 410 411 // Marks results 412 const std::vector<double>& X = solver.solution(); 413 414 std::vector<Face_descriptor> to_delete; 415 std::size_t f_idx(0); 416 for(auto f : target_mesh.faces()) { 417 if (static_cast<int>(std::round(X[f_idx])) == 0) 418 to_delete.push_back(f); 419 ++f_idx; 420 } 421 422 for (std::size_t i = 0; i < to_delete.size(); ++i) { 423 Face_descriptor f = to_delete[i]; 424 Halfedge_descriptor h = target_mesh.halfedge(f); 425 Euler::remove_face(h, target_mesh); 426 } 427 428 // Marks the sharp edges 429 typename Polygon_mesh::template Property_map<Edge_descriptor, bool> edge_is_sharp = 430 target_mesh.template add_property_map<Edge_descriptor, bool>("e:sharp_edges").first; 431 for (auto e : target_mesh.edges()) 432 edge_is_sharp[e] = false; 433 434 for (std::size_t i = 0; i < adjacency.size(); ++i) { 435 const Intersection& fan = adjacency[i]; 436 if (fan.size() != 4) 437 continue; 438 439 std::size_t idx_sharp_var = edge_sharp_status[&fan]; 440 if (static_cast<int>(X[idx_sharp_var]) == 1) { 441 for (std::size_t j = 0; j < fan.size(); ++j) { 442 Halfedge_descriptor h = fan[j]; 443 Face_descriptor f = target_mesh.face(h); 444 if (f != Polygon_mesh::null_face()) { // some faces may be deleted 445 std::size_t fid = face_indices[f]; 446 if (static_cast<int>(std::round(X[fid])) == 1) { 447 Edge_descriptor e = target_mesh.edge(h); 448 edge_is_sharp[e] = true; 449 break; 450 } 451 } 452 } 453 } 454 } 455 456 // Converts from internal data structure to the required `PolygonMesh`. 457 output_mesh.clear(); // make sure it is empty. 458 CGAL::copy_face_graph(target_mesh, output_mesh); 459 } 460 else { 461 error_message_ = "solving the binary program failed"; 462 return false; 463 } 464 465 return true; 466 } 467 468 } //namespace CGAL 469 470 #endif // CGAL_POLYGONAL_SURFACE_RECONSTRUCTION_H 471
