diff --git a/modules/odm_25dmeshing/src/Odm25dMeshing.cpp b/modules/odm_25dmeshing/src/Odm25dMeshing.cpp index 3e16c037..23651bb8 100644 --- a/modules/odm_25dmeshing/src/Odm25dMeshing.cpp +++ b/modules/odm_25dmeshing/src/Odm25dMeshing.cpp @@ -1,25 +1,14 @@ #include "Odm25dMeshing.hpp" //We define a vertex_base with info. The "info" (size_t) allow us to keep track of the original point index. -//typedef CGAL::Triangulation_vertex_base_with_info_2 Vb; -//typedef CGAL::Triangulation_data_structure_2 Tds; -//typedef CGAL::Delaunay_triangulation_2 DT; -//typedef DT::Point cgalPoint; -//typedef DT::Vertex_circulator Vertex_circulator; +typedef CGAL::Triangulation_vertex_base_with_info_2 Vb; +typedef CGAL::Triangulation_data_structure_2 Tds; +typedef CGAL::Delaunay_triangulation_2 DT; +typedef DT::Point cgalPoint; +typedef DT::Vertex_circulator Vertex_circulator; -typedef CGAL::Projection_traits_xy_3 Gt; -typedef CGAL::Triangulation_vertex_base_2 Vb; -typedef CGAL::Delaunay_mesh_face_base_2 Fb; -typedef CGAL::Triangulation_data_structure_2 Tds; - -typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; - -typedef CGAL::Delaunay_mesh_size_criteria_2 Criteria; -typedef CGAL::Delaunay_mesher_2 Mesher; - -typedef CDT::Vertex_handle Vertex_handle; -typedef CDT::Vertex_circulator Vertex_circulator; -typedef CDT::Point cgalPoint; +typedef CGAL::First_of_pair_property_map Point_map; +typedef CGAL::Second_of_pair_property_map Normal_map; // Concurrency #ifdef CGAL_LINKED_WITH_TBB @@ -41,7 +30,7 @@ int Odm25dMeshing::run(int argc, char **argv) { parseArguments(argc, argv); -// loadPointCloud(); + loadPointCloud(); buildMesh(); @@ -83,15 +72,6 @@ void Odm25dMeshing::parseArguments(int argc, char **argv) { if (ss.bad()) throw Odm25dMeshingException("Argument '" + argument + "' has a bad value (wrong type)."); maxVertexCount = std::max(maxVertexCount, 0); log << "Vertex count was manually set to: " << maxVertexCount << "\n"; - } else if (argument == "-outliersRemovalPercentage" && argIndex < argc) { - ++argIndex; - if (argIndex >= argc) throw Odm25dMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided."); - std::stringstream ss(argv[argIndex]); - ss >> outliersRemovalPercentage; - if (ss.bad()) throw Odm25dMeshingException("Argument '" + argument + "' has a bad value (wrong type)."); - - outliersRemovalPercentage = std::min(99.99, std::max(outliersRemovalPercentage, 0)); - log << "Outliers removal was manually set to: " << outliersRemovalPercentage << "\n"; } else if (argument == "-wlopIterations" && argIndex < argc) { ++argIndex; if (argIndex >= argc) throw Odm25dMeshingException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided."); @@ -167,90 +147,53 @@ void Odm25dMeshing::loadPointCloud(){ } void Odm25dMeshing::buildMesh(){ -// const unsigned int NEIGHBORS = 24; -// -// if (outliersRemovalPercentage > 0){ -// size_t pointCountBeforeOutRemoval = points.size(); -// log << "Removing outliers\n"; -// -// points.erase(CGAL::remove_outliers(points.begin(), points.end(), -// CGAL::First_of_pair_property_map(), NEIGHBORS, outliersRemovalPercentage), -// points.end()); -// std::vector(points).swap(points); -// -// log << "Removed " << (pointCountBeforeOutRemoval - points.size()) << " points\n"; -// } -// -// // Filter points by normal orientation -// log << "Applying normal filtering\n"; -// -// double COSINE_LIMIT = std::cos(M_PI / 4); // ~ 45 degrees -// -// Vector3 up(0, 0, 1); -// if (flipFaces) up = Vector3(0, 0, -1); -// -// size_t pointCountBeforeNormalFilter = points.size(); -// -// points.erase(std::remove_if(points.begin(), points.end(),[&](Pwn pointWithNormal){ -// double cosine = pointWithNormal.second * up / CGAL::sqrt(pointWithNormal.second*pointWithNormal.second) / CGAL::sqrt(up * up); -// return cosine <= COSINE_LIMIT; -// }), points.end()); -// + const unsigned int NEIGHBORS = 16; + size_t pointCount = points.size(); -// -// log << "Removed " << (pointCountBeforeNormalFilter - pointCount) << " points\n"; -// -// const double RETAIN_PERCENTAGE = std::min(((100. * (double)maxVertexCount) / (double)pointCount), 80.); // percentage of points to retain. - std::vector simplifiedPoints; -// -// log << "Computing average spacing... "; -// -// FT avgSpacing = CGAL::compute_average_spacing( -// points.begin(), -// points.end(), -// CGAL::First_of_pair_property_map(), -// NEIGHBORS); -// -// log << avgSpacing << "\n"; -// -// log << "Performing weighted locally optimal projection simplification and regularization (retain: " << RETAIN_PERCENTAGE << "%, iterate: " << wlopIterations << ")" << "\n"; -// -// CGAL::wlop_simplify_and_regularize_point_set( -// points.begin(), -// points.end(), -// std::back_inserter(simplifiedPoints), -// CGAL::First_of_pair_property_map(), -// RETAIN_PERCENTAGE, -// 8 * avgSpacing, -// wlopIterations, -// false); // require_uniform_sampling -// for (size_t c = 0; c < pointCount; c++){ -// simplifiedPoints.push_back(points[c].first); -// } - -// std::ofstream f("out.bin"); -// f << simplifiedPoints.size() * 3 << " "; -// -// for (size_t i = 0; i < simplifiedPoints.size(); i++){ -// f << simplifiedPoints[i].x() << " "; -// f << simplifiedPoints[i].y() << " "; -// f << simplifiedPoints[i].z() << " "; -// } -// f.close(); -// exit(1); - - std::ifstream fin("out.bin"); - pointCount = 0; - fin >> pointCount; - for (size_t i = 0; i < pointCount / 3; i++){ - float x, y, z; - fin >> x; - fin >> y; - fin >> z; - simplifiedPoints.push_back(Point3(x, y, z)); + log << "Performing bilateral smoothing... "; + for (int i = 0; i < 3; i++){ + double err = CGAL::bilateral_smooth_point_set( + points.begin(), + points.end(), + Point_map(), + Normal_map(), + NEIGHBORS, + 45.0, + Kernel()); + log << err << " "; + if (err < 0.001){ + log << "stopping early\n"; + }else if (i >= 2){ + log << "iteration limit reached\n"; + } } + const double RETAIN_PERCENTAGE = std::min(((100. * (double)maxVertexCount) / (double)pointCount), 80.); // percentage of points to retain. + std::vector simplifiedPoints; + + log << "Computing average spacing... "; + + FT avgSpacing = CGAL::compute_average_spacing( + points.begin(), + points.end(), + Point_map(), + NEIGHBORS); + + log << avgSpacing << "\n"; + + log << "Performing weighted locally optimal projection simplification and regularization (retain: " << RETAIN_PERCENTAGE << "%, iterate: " << wlopIterations << ")" << "\n"; + + CGAL::wlop_simplify_and_regularize_point_set( + points.begin(), + points.end(), + std::back_inserter(simplifiedPoints), + Point_map(), + RETAIN_PERCENTAGE, + 8 * avgSpacing, + wlopIterations, + true); + pointCount = simplifiedPoints.size(); if (pointCount < 3){ @@ -259,21 +202,74 @@ void Odm25dMeshing::buildMesh(){ log << "Vertex count is " << pointCount << "\n"; -// std::vector< std::pair > pts; -// try{ -// pts.reserve(pointCount); -// } catch (const std::bad_alloc&){ -// throw Odm25dMeshingException("Not enough memory"); -// } + log << "Z-occlusion filtering\n"; + + size_t pointCountBeforeOcclusionFiltering = pointCount; + + const double gridStep = 0.05; + Kernel::Iso_cuboid_3 bbox = CGAL::bounding_box(simplifiedPoints.begin(), simplifiedPoints.end()); + Vector3 boxDiag = bbox.max() - bbox.min(); + + int gridWidth = 1 + static_cast(boxDiag.x() / gridStep + 0.5); + int gridHeight = 1 + static_cast(boxDiag.y() / gridStep + 0.5); + + std::unordered_map grid; + + for (size_t c = 0; c < pointCount; c++){ + const Point3 &p = simplifiedPoints[c]; + Vector3 relativePos = p - bbox.min(); + int i = static_cast((relativePos.x() / gridStep + 0.5)); + int j = static_cast((relativePos.y() / gridStep + 0.5)); + + if ((i >= 0 && i < gridWidth) && (j >= 0 && j < gridHeight)){ + int key = i * gridWidth + j; + + if (grid.find(key) == grid.end()){ + grid[key] = p; + }else if (p.z() > grid[key].z()){ + grid[key] = p; + } // -// for (size_t i = 0; i < pointCount; ++i){ -// pts.push_back(std::make_pair(cgalPoint(simplifiedPoints[i].x(), simplifiedPoints[i].y()), i)); +// if (grid[i][j] == NULL){ +// grid[i][j] = static_cast(&p); +// }else if (p.z() > grid[i][j]->z()){ +// grid[i][j] = static_cast(&p); +// } + } + } + + std::vector gridPoints; + for ( auto it = grid.begin(); it != grid.end(); ++it ){ + gridPoints.push_back(it->second); + } + +// for (int i = 0; i < gridWidth; i++){ +// for (int j = 0; j < gridHeight; j++){ +// if (grid[i][j] != NULL){ +// gridPoints.push_back(*grid[i][j]); +// } +// } // } + pointCount = gridPoints.size(); + log << "Removed " << (pointCountBeforeOcclusionFiltering - pointCount) << " points\n"; + + + std::vector< std::pair > pts; + try{ + pts.reserve(pointCount); + } catch (const std::bad_alloc&){ + throw Odm25dMeshingException("Not enough memory"); + } + + for (size_t i = 0; i < pointCount; ++i){ + pts.push_back(std::make_pair(cgalPoint(gridPoints[i].x(), gridPoints[i].y()), i)); + } + log << "Computing delaunay triangulation\n"; //The delaunay triangulation is built according to the 2D point cloud -// DT dt(pts.begin(), pts.end()); + DT dt(pts.begin(), pts.end()); unsigned int numberOfTriangles = static_cast(dt.number_of_faces()); unsigned int triIndexes = dt.number_of_faces()*3; @@ -295,9 +291,9 @@ void Odm25dMeshing::buildMesh(){ for (size_t i = 0; i < pointCount; ++i){ - vertices.push_back(simplifiedPoints[i].x()); - vertices.push_back(simplifiedPoints[i].y()); - vertices.push_back(simplifiedPoints[i].z()); + vertices.push_back(gridPoints[i].x()); + vertices.push_back(gridPoints[i].y()); + vertices.push_back(gridPoints[i].z()); } for (DT::Face_iterator face = dt.faces_begin(); face != dt.faces_end(); ++face) { @@ -314,12 +310,11 @@ void Odm25dMeshing::buildMesh(){ log << "Removing spikes\n"; - const float THRESHOLD = 0.1; + const float THRESHOLD = 0.05; std::vector heights; unsigned int spikesRemoved = 0; - for (DT::Vertex_iterator vertex = dt.vertices_begin(); vertex != dt.vertices_end(); ++vertex){ // Check if the height between this vertex and its // incident vertices is greater than THRESHOLD @@ -387,7 +382,6 @@ void Odm25dMeshing::printHelp() { << " -outputFile where the output PLY 2.5D mesh should be saved (default: " << outputFile << ")\n" << " -logFile log file path (default: " << logFilePath << ")\n" << " -verbose whether to print verbose output (default: " << (printInCoutPop ? "true" : "false") << ")\n" - << " -outliersRemovalPercentage <0 - 99.99> percentage of outliers to remove. Set to 0 to disable. (default: " << outliersRemovalPercentage << ")\n" << " -maxVertexCount <0 - N> Maximum number of vertices in the output mesh. The mesh might have fewer vertices, but will not exceed this limit. (default: " << maxVertexCount << ")\n" << " -wlopIterations <1 - 1000> Iterations of the Weighted Locally Optimal Projection (WLOP) simplification algorithm. Higher values take longer but produce a smoother mesh. (default: " << wlopIterations << ")\n" diff --git a/modules/odm_25dmeshing/src/Odm25dMeshing.hpp b/modules/odm_25dmeshing/src/Odm25dMeshing.hpp index fd6d9ab8..c57c2356 100644 --- a/modules/odm_25dmeshing/src/Odm25dMeshing.hpp +++ b/modules/odm_25dmeshing/src/Odm25dMeshing.hpp @@ -3,17 +3,16 @@ // STL #include #include +#include #include #include #include #include #include -#include -#include #include -#include -#include +#include +#include #include "Logger.hpp" #include "PlyInterpreter.hpp" @@ -65,7 +64,6 @@ private: std::string outputFile = "odm_25dmesh.ply"; std::string logFilePath = "odm_25dmeshing_log.txt"; unsigned int maxVertexCount = 100000; - double outliersRemovalPercentage = 2; unsigned int wlopIterations = 35; std::vector points; bool flipFaces = false; diff --git a/opendm/config.py b/opendm/config.py index b1683782..9cd7653c 100644 --- a/opendm/config.py +++ b/opendm/config.py @@ -262,14 +262,6 @@ def config(): 'Applies to 2.5D mesh only. ' 'Default: %(default)s')) - parser.add_argument('--mesh-remove-outliers', - metavar='', - default=2, - type=float, - help=('Percentage of outliers to remove from the point set. Set to 0 to disable. ' - 'Applies to 2.5D mesh only. ' - 'Default: %(default)s')) - parser.add_argument('--mesh-wlop-iterations', metavar='', default=70, diff --git a/scripts/odm_app.py b/scripts/odm_app.py index fabe2463..6c6c481a 100644 --- a/scripts/odm_app.py +++ b/scripts/odm_app.py @@ -59,7 +59,6 @@ class ODMApp(ecto.BlackBox): samples=p.args.mesh_samples, solver=p.args.mesh_solver_divide, max_vertex_25d=getattr(p.args, '25d_mesh_size'), - remove_outliers=p.args.mesh_remove_outliers, wlop_iterations=p.args.mesh_wlop_iterations, verbose=p.args.verbose), 'texturing': ODMMvsTexCell(data_term=p.args.texturing_data_term, diff --git a/scripts/odm_meshing.py b/scripts/odm_meshing.py index e02fb4c7..32e91e9a 100644 --- a/scripts/odm_meshing.py +++ b/scripts/odm_meshing.py @@ -21,9 +21,6 @@ class ODMeshingCell(ecto.Cell): 'times slightly but helps reduce memory usage.', 9) params.declare("max_vertex_25d", 'The maximum vertex count of the 2.5D output mesh.', 50000) - params.declare("remove_outliers", 'Percentage of outliers to remove from the point set. Set to 0 to disable. ' - 'Applies to 2.5D mesh only.', 2) - params.declare("wlop_iterations", 'Iterations of the Weighted Locally Optimal Projection (WLOP) simplification algorithm. ' 'Higher values take longer but produce a smoother mesh. ' 'Applies to 2.5D mesh only. ', 70) @@ -98,14 +95,13 @@ class ODMeshingCell(ecto.Cell): 'log': tree.odm_25dmeshing_log, 'verbose': verbose, 'max_vertex': self.params.max_vertex_25d, - 'remove_outliers': self.params.remove_outliers, 'wlop_iterations': self.params.wlop_iterations } # run 2.5D meshing binary system.run('{bin}/odm_25dmeshing -inputFile {infile} ' '-outputFile {outfile} -logFile {log} ' - '-maxVertexCount {max_vertex} -outliersRemovalPercentage {remove_outliers} ' + '-maxVertexCount {max_vertex} ' '-wlopIterations {wlop_iterations} {verbose}'.format(**kwargs)) else: log.ODM_WARNING('Found a valid ODM 2.5D Mesh file in: %s' %