cgal点云重建白膜_cejutue的博客-爱代码爱编程

基本效果如下:

主要代码如下:

#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库