基本效果如下:
主要代码如下:
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_xyz_points.h>
#include <CGAL/IO/Writer_OFF.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/property_map.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
#include <CGAL/Polygonal_surface_reconstruction.h>
#include <CGAL/remove_outliers.h>
#include <CGAL/grid_simplify_point_set.h>
#include <CGAL/compute_average_spacing.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/jet_estimate_normals.h>
#include <CGAL/pca_estimate_normals.h>
#include <CGAL/mst_orient_normals.h>
#include <CGAL/radial_orient_normals.h>
#include <CGAL/vcm_estimate_normals.h>
#define CGAL_USE_GLPK
#ifdef CGAL_USE_SCIP // defined (or not) by CMake scripts, do not define by hand
#include <CGAL/SCIP_mixed_integer_program_traits.h>
typedef CGAL::SCIP_mixed_integer_program_traits<double> MIP_Solver;
#elif defined(CGAL_USE_GLPK) // defined (or not) by CMake scripts, do not define by hand
#include <CGAL/GLPK_mixed_integer_program_traits.h>
typedef CGAL::GLPK_mixed_integer_program_traits<double> MIP_Solver;
#endif
#if defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)
#include <CGAL/Timer.h>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;
// Point with normal, and plane index
typedef boost::tuple<Point, Vector, int> PNI;
typedef std::vector<PNI> Point_vector;
//using Point_vector = CGAL::Point_set_3<PNI>;
typedef CGAL::Nth_of_tuple_property_map<0, PNI> Point_map;
typedef CGAL::Nth_of_tuple_property_map<1, PNI> Normal_map;
typedef CGAL::Nth_of_tuple_property_map<2, PNI> Plane_index_map;
typedef CGAL::Shape_detection::Efficient_RANSAC_traits<Kernel, Point_vector, Point_map, Normal_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection::Plane<Traits> Plane;
typedef CGAL::Shape_detection::Sphere<Traits> Sphere;
typedef CGAL::Shape_detection::Cone<Traits> Cone;
typedef CGAL::Shape_detection::Torus<Traits> Torus;
typedef CGAL::Shape_detection::Cylinder<Traits> Cylinder;
typedef CGAL::Shape_detection::Point_to_shape_index_map<Traits> Point_to_shape_index_map;
typedef CGAL::Polygonal_surface_reconstruction<Kernel> Polygonal_surface_reconstruction;
typedef CGAL::Surface_mesh<Point> Surface_mesh;
/*
* This example first extracts planes from the input point cloud
* (using RANSAC with default parameters) and then reconstructs
* the surface model from the planes.
*/
/*
*/
// Boost includes.
#include <boost/function_output_iterator.hpp>
// CGAL includes.
#include <CGAL/Timer.h>
#include <CGAL/Random.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Point_set_3/IO.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
#include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h>
#include <CGAL/jet_smooth_point_set.h>
class Index_map {
public:
using key_type = std::size_t;
using value_type = int;
using reference = value_type;
using category = boost::readable_property_map_tag;
Index_map() { }
template<typename PointRange>
Index_map(
const PointRange& points,
const std::vector< std::vector<std::size_t> >& regions) :
m_indices(new std::vector<int>(points.size(), -1)) {
for (std::size_t i = 0; i < regions.size(); ++i)
for (const std::size_t idx : regions[i])
(*m_indices)[idx] = static_cast<int>(i);
(*m_indices)[points.size() - 1] = regions.size();
(*m_indices)[points.size() - 2] = regions.size();
}
inline friend value_type get(
const Index_map& index_map,
const key_type key) {
const auto& indices = *(index_map.m_indices);
return indices[key];
}
void bujiu()
{
}
public:
std::shared_ptr< std::vector<int> > m_indices;
};
int region_grow(std::string & filename)
{
CGAL::Timer t;
// Type declarations.
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = typename Kernel::FT;
using Point_3 = typename Kernel::Point_3;
using Vector_3 = typename Kernel::Vector_3;
using Input_range = CGAL::Point_set_3<Point_3>;
using Point_map = typename Input_range::Point_map;
using Normal_map = typename Input_range::Vector_map;
using Neighbor_query = CGAL::Shape_detection::Point_set::K_neighbor_query<Kernel, Input_range, Point_map>;
using Region_type = CGAL::Shape_detection::Point_set::Least_squares_plane_fit_region<Kernel, Input_range, Point_map, Normal_map>;
using Region_growing = CGAL::Shape_detection::Region_growing<Input_range, Neighbor_query, Region_type>;
using Indices = std::vector<std::size_t>;
using Output_range = CGAL::Point_set_3<Point_3>;
using Points_3 = std::vector<Point_3>;
using Traits = CGAL::Shape_detection::Efficient_RANSAC_traits
<Kernel, Input_range, Input_range::Point_map, Input_range::Vector_map>;
using Efficient_ransac = CGAL::Shape_detection::Efficient_RANSAC<Traits>;
using Plane = CGAL::Shape_detection::Plane<Traits>;
using Sphere = CGAL::Shape_detection::Sphere<Traits>;
using Cone = CGAL::Shape_detection::Cone<Traits>;
using Torus = CGAL::Shape_detection::Torus<Traits>;
using Cylinder = CGAL::Shape_detection::Cylinder<Traits>;
using Point_to_shape_index_map = CGAL::Shape_detection::Point_to_shape_index_map<Traits>;
using Polygonal_surface_reconstruction = CGAL::Polygonal_surface_reconstruction<Kernel>;
using Surface_mesh = CGAL::Surface_mesh<Point>;
// Concurrency
using Concurrency_tag = CGAL::Parallel_if_available_tag;
// Define an insert iterator.
struct Insert_point_colored_by_region_index {
using argument_type = Indices;
using result_type = void;
using Color_map =
typename Output_range:: template Property_map<unsigned char>;
const Input_range& m_input_range;
const Point_map m_point_map;
Output_range& m_output_range;
std::size_t& m_number_of_regions;
Color_map m_red, m_green, m_blue;
Insert_point_colored_by_region_index(
const Input_range& input_range,
const Point_map point_map,
Output_range& output_range,
std::size_t& number_of_regions) :
m_input_range(input_range),
m_point_map(point_map),
m_output_range(output_range),
m_number_of_regions(number_of_regions) {
m_red =
m_output_range.template add_property_map<unsigned char>("red", 0).first;
m_green =
m_output_range.template add_property_map<unsigned char>("green", 0).first;
m_blue =
m_output_range.template add_property_map<unsigned char>("blue", 0).first;
}
result_type operator()(const argument_type& region) {
CGAL::Random rand(static_cast<unsigned int>(m_number_of_regions));
const unsigned char r =
static_cast<unsigned char>(64 + rand.get_int(0, 192));
const unsigned char g =
static_cast<unsigned char>(64 + rand.get_int(0, 192));
const unsigned char b =
static_cast<unsigned char>(64 + rand.get_int(0, 192));
for (const std::size_t index : region) {
const auto& key = *(m_input_range.begin() + index);
const Point_3& point = get(m_point_map, key);
const auto it = m_output_range.insert(point);
m_red[*it] = r;
m_green[*it] = g;
m_blue[*it] = b;
}
++m_number_of_regions;
}
}; // Insert_point_colored_by_region_index
// Load xyz data either from a local folder or a user-provided file.
//std::ifstream in("data/Tile_+015_+029-twohouser.xyz");
std::string pathname = filename;//"D:\\testdata\\xyz\\Tile_+014_+029_1.xyz";
std::ifstream in(pathname);
// Reading input in XYZ format
Input_range point_set;
if (!in || !CGAL::read_xyz_point_set(in, point_set))
{
std::cerr << "Can't read input file " << std::endl;
return 0;
}
//Simplify point set
//auto points_2_plane = [](Input_range& Input,std::vector< std::vector<std::size_t> > regions) {
//
//};
double spacing = CGAL::compute_average_spacing<CGAL::Sequential_tag>(point_set, 70);
auto gsim_it = CGAL::grid_simplify_point_set(point_set, spacing);
point_set.remove(gsim_it, point_set.end());
CGAL::jet_smooth_point_set<CGAL::Sequential_tag>(point_set, 24);
point_set.collect_garbage();
auto box = CGAL::bounding_box(point_set.points().begin(), point_set.points().end());
Point pminxy(box.xmin(), box.ymin(), box.zmin());
Point pminz(box.xmax(), box.ymax(), box.zmin());
point_set.insert(pminxy);
point_set.insert(pminz);
//if (point_set.has_normal_map())
//{
// for (Input_range::iterator it = point_set.begin(); it != point_set.end(); ++it)
// {
// Vector n = point_set.normal(*it);
// n = -n / std::sqrt(n * n);
// point_set.normal(*it) = n;
// }
//}
//else
//{
if (point_set.has_normal_map())
point_set.remove_normal_map();
point_set.add_normal_map();
CGAL::jet_estimate_normals<CGAL::Sequential_tag>(point_set, 8, point_set.parameters().degree_fitting(2)); // additional named parameter specific to jet_estimate_normals
//}
if (!point_set.has_normal_map())
{
std::cerr << " Failed: has_normal_map" << std::endl;
return 0;
}
std::cout << "* loaded " << point_set.size() << " points with normals" << std::endl;
// Default parameter values for the data file point_set_3.xyz.
const std::size_t k = 16;
const FT max_distance_to_plane = FT(spacing * 1.21);
const FT max_accepted_angle = FT(35);
const std::size_t min_region_size = 50;
// Create instances of the classes Neighbor_query and Region_type.
Neighbor_query neighbor_query(
point_set,
k,
point_set.point_map());
Region_type region_type(
point_set,
max_distance_to_plane, max_accepted_angle, min_region_size,
point_set.point_map(), point_set.normal_map());
// Create an instance of the region growing class.
Region_growing region_growing(
point_set, neighbor_query, region_type);
// Run the algorithm.
Output_range output_range;
CGAL::Timer timer;
std::vector< std::vector<std::size_t> > regions;
timer.start();
region_growing.detect(std::back_inserter(regions));
std::size_t number_of_regions = 0;
Insert_point_colored_by_region_index inserter(
point_set, point_set.point_map(),
output_range, number_of_regions);
region_growing.detect(
boost::make_function_output_iterator(inserter));
timer.stop();
// Print the number of found regions.
std::cout << "* " << regions.size() << " regions have been found in " << timer.time() << " seconds"
<< std::endl;
std::string tmpply = pathname;
tmpply = tmpply.replace(tmpply.find(".xyz"), 4, ".ply");
const std::string path = tmpply;//"data/Tile_+015_+029-twohouser3333.ply";
std::ofstream out(path);
out << output_range;
std::cout << "* found regions are saved in " << path << std::endl;
out.close();
Index_map index_map(point_set, regions);
Polygonal_surface_reconstruction algo(
point_set,
point_set.point_map(),
point_set.normal_map(),
index_map);
Surface_mesh model;
std::cout << "Reconstructing...";
if (!algo.reconstruct<MIP_Solver>(model, 0.43, 0.27, 0.3)) {
std::cerr << " Failed: " << algo.error_message() << std::endl;
return 0;
}
pathname = pathname.replace(pathname.find(".xyz"), 4, ".off");
const std::string& output_file(pathname);
std::ofstream output_stream(output_file.c_str());
if (output_stream && CGAL::write_off(output_stream, model)) {
// flush the buffer
output_stream << std::flush;
std::cout << " Done. Saved to " << output_file << ". Time: " << t.time() << " sec." << std::endl;
}
else {
std::cerr << " Failed saving file." << std::endl;
return 0;
}
return 1;
}
此方案需要先将房子一栋栋分开, 目前没研究, 房子是手动分开, 一栋栋处理的, CGAL使用版本5.1.5,boost6.18,gplk库