diff --git a/modules/odm_25dmeshing/src/Odm25dMeshing.cpp b/modules/odm_25dmeshing/src/Odm25dMeshing.cpp index 23651bb8..3459bce4 100644 --- a/modules/odm_25dmeshing/src/Odm25dMeshing.cpp +++ b/modules/odm_25dmeshing/src/Odm25dMeshing.cpp @@ -7,8 +7,8 @@ typedef CGAL::Delaunay_triangulation_2 DT; typedef DT::Point cgalPoint; typedef DT::Vertex_circulator Vertex_circulator; -typedef CGAL::First_of_pair_property_map Point_map; -typedef CGAL::Second_of_pair_property_map Normal_map; +//typedef CGAL::First_of_pair_property_map Point_map; +//typedef CGAL::Second_of_pair_property_map Normal_map; // Concurrency #ifdef CGAL_LINKED_WITH_TBB @@ -147,48 +147,101 @@ void Odm25dMeshing::loadPointCloud(){ } void Odm25dMeshing::buildMesh(){ - const unsigned int NEIGHBORS = 16; + const unsigned int NEIGHBORS = 24; size_t pointCount = points.size(); - 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 << "Grid Z sampling\n"; + + size_t pointCountBeforeGridSampling = pointCount; + + double gridStep = avgSpacing; + Kernel::Iso_cuboid_3 bbox = CGAL::bounding_box(points.begin(), points.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); + + #define KEY(i, j) (i * gridWidth + j) + + std::unordered_map grid; + + for (size_t c = 0; c < pointCount; c++){ + const Point3 &p = points[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 = KEY(i, j); + + if (grid.find(key) == grid.end()){ + grid[key] = p; + }else if ((!flipFaces && p.z() > grid[key].z()) || (flipFaces && p.z() < grid[key].z())){ + grid[key] = p; + } + } + } + + std::vector gridPoints; + for ( auto it = grid.begin(); it != grid.end(); ++it ){ + gridPoints.push_back(it->second); + } + + pointCount = gridPoints.size(); + log << "Sampled " << (pointCountBeforeGridSampling - pointCount) << " points\n"; + + log << "Hole filling\n"; + + std::unordered_map::iterator it; + for (int i = 1; i < gridWidth - 1; i++){ + for (int j = 1; j < gridHeight - 1; j++){ + if (grid.find(KEY(i, j)) == grid.end()){ + int neighbors = 0; + FT avgZ = 0; + + for (int k = i - 1; k < i + 1; k++){ + for (int l = j - 1; l < j + 1; l++){ + // Search for immediate neighbors + if ((it = grid.find(KEY(k, l))) != grid.end()){ + avgZ += it->second.z(); + neighbors++; + } + } + } + + // Found? + if (neighbors > 0){ + gridPoints.push_back(Point3( + (i - 0.5) * gridStep + bbox.min().x(), + (j - 0.5) * gridStep + bbox.min().y(), + avgZ / static_cast(neighbors))); + } + } + } + } + + log << "Filled grid with " << (gridPoints.size() - pointCount) << " points\n"; + pointCount = gridPoints.size(); + + const double RETAIN_PERCENTAGE = std::min(100., 100. * static_cast(maxVertexCount) / static_cast(pointCount)); // percentage of points to retain. + std::vector simplifiedPoints; + 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(), + gridPoints.begin(), + gridPoints.end(), std::back_inserter(simplifiedPoints), - Point_map(), RETAIN_PERCENTAGE, 8 * avgSpacing, wlopIterations, @@ -202,59 +255,6 @@ void Odm25dMeshing::buildMesh(){ log << "Vertex count is " << pointCount << "\n"; - 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; - } -// -// 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); @@ -263,7 +263,7 @@ void Odm25dMeshing::buildMesh(){ } for (size_t i = 0; i < pointCount; ++i){ - pts.push_back(std::make_pair(cgalPoint(gridPoints[i].x(), gridPoints[i].y()), i)); + pts.push_back(std::make_pair(cgalPoint(simplifiedPoints[i].x(), simplifiedPoints[i].y()), i)); } log << "Computing delaunay triangulation\n"; @@ -291,9 +291,9 @@ void Odm25dMeshing::buildMesh(){ for (size_t i = 0; i < pointCount; ++i){ - vertices.push_back(gridPoints[i].x()); - vertices.push_back(gridPoints[i].y()); - vertices.push_back(gridPoints[i].z()); + vertices.push_back(simplifiedPoints[i].x()); + vertices.push_back(simplifiedPoints[i].y()); + vertices.push_back(simplifiedPoints[i].z()); } for (DT::Face_iterator face = dt.faces_begin(); face != dt.faces_end(); ++face) { diff --git a/modules/odm_25dmeshing/src/Odm25dMeshing.hpp b/modules/odm_25dmeshing/src/Odm25dMeshing.hpp index c57c2356..371e03b6 100644 --- a/modules/odm_25dmeshing/src/Odm25dMeshing.hpp +++ b/modules/odm_25dmeshing/src/Odm25dMeshing.hpp @@ -65,7 +65,7 @@ private: std::string logFilePath = "odm_25dmeshing_log.txt"; unsigned int maxVertexCount = 100000; unsigned int wlopIterations = 35; - std::vector points; + std::vector points; bool flipFaces = false; }; diff --git a/modules/odm_25dmeshing/src/PlyInterpreter.cpp b/modules/odm_25dmeshing/src/PlyInterpreter.cpp index c8e7fc55..f8a07488 100644 --- a/modules/odm_25dmeshing/src/PlyInterpreter.cpp +++ b/modules/odm_25dmeshing/src/PlyInterpreter.cpp @@ -19,11 +19,11 @@ void PlyInterpreter::process_line(CGAL::Ply_reader& reader) { reader.assign (y, "y"); reader.assign (z, "z"); reader.assign (nx, "nx"); - reader.assign (ny, "ny"); - reader.assign (nz, "nz"); + reader.assign (ny, "ny"); + reader.assign (nz, "nz"); Point3 p(x, y, z); - Vector3 n(nx, ny, nz); +// Vector3 n(nx, ny, nz); if (nz >= 0 && zNormalsDirectionCount < std::numeric_limits::max()){ zNormalsDirectionCount++; @@ -31,7 +31,8 @@ void PlyInterpreter::process_line(CGAL::Ply_reader& reader) { zNormalsDirectionCount--; } - points.push_back(std::make_pair(p, n)); +// points.push_back(std::make_pair(p, n)); + points.push_back(p); } bool PlyInterpreter::flip_faces(){ diff --git a/modules/odm_25dmeshing/src/PlyInterpreter.hpp b/modules/odm_25dmeshing/src/PlyInterpreter.hpp index f7b78478..3d0eeaab 100644 --- a/modules/odm_25dmeshing/src/PlyInterpreter.hpp +++ b/modules/odm_25dmeshing/src/PlyInterpreter.hpp @@ -16,14 +16,14 @@ typedef Kernel::Point_3 Point3; typedef Kernel::Vector_3 Vector3; // points, normals -typedef std::pair Pwn; +//typedef std::pair Pwn; class PlyInterpreter { - std::vector& points; + std::vector& points; long zNormalsDirectionCount; public: - PlyInterpreter (std::vector& points) + PlyInterpreter (std::vector& points) : points (points), zNormalsDirectionCount(0) { } bool is_applicable (CGAL::Ply_reader& reader); diff --git a/opendm/config.py b/opendm/config.py index 2c3b9afe..1fe060d5 100644 --- a/opendm/config.py +++ b/opendm/config.py @@ -254,7 +254,7 @@ def config(): parser.add_argument('--mesh-wlop-iterations', metavar='', - default=70, + default=35, type=int, help=('Iterations of the Weighted Locally Optimal Projection (WLOP) simplification algorithm. ' 'Higher values take longer but produce a smoother mesh. '