给定点，计算法线并在 CGAL 中执行多边形表面重建

问题描述

我有 .xyz 文件，我想在 Cgal 中使用“多边形表面重建”重新创建表面。重建表面需要法线，所以我想首先在 Cgal 中计算这些，然后使用这些点及其法线作为多边形表面重建的输入。

但是，两者的输出和输入类型不同：

正常估计的输出是一个vectorpairs列表，见https://doc.cgal.org/latest/Point_set_processing_3/index.html#Point_set_processing_3NormalOrientation：12.1.1示例；而多边形表面重建的输入是一个 point_vector ((https://doc.cgal.org/latest/Polygonal_surface_reconstruction/index.html).

我尝试创建一个转换数据格式的“convert_function”，但我不是 C++ 程序员，因此我很难转换。如果有人可以帮助我解决这个问题，那就太好了，因为我已经坚持了几个小时。我的目标是让这个过程尽可能高效，因此非常感谢关于如何解决这个问题的更好的想法。

数据类型：

->正常计算：

typedef std::pair<Point,Vector> PointVectorPair;

->表面重建：

typedef boost::tuple<Point,Vector,int> PNI;
typedef std::vector<PNI> Point_vector;
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;

我的反转功能不起作用，因为语法不正确：

std::vector<PNI> convert_function(const std::list<PointVectorPair>& list)
{
    std::vector<PNI> out;//a temporary object to store the output

    for (std::size_t i = 0; i < list.size(); ++i) {
    {   
        out[i].get<0>() = list[i].get<0>();
        out[i].get<1>() = list[i].get<1>();

    }
    return out;
}

我再进一步使用以下内容：

std::vector<PNI> points2; // store points
points2 = convert_function(points);

现在完整的代码如下：

#include <Cgal/Exact_predicates_inexact_constructions_kernel.h>
#include <Cgal/IO/read_xyz_points.h>
#include <Cgal/IO/Writer_OFF.h>
#include <Cgal/property_map.h>
#include <Cgal/Surface_mesh.h>
#include <Cgal/Shape_detection/Efficient_RANSAC.h>
#include <Cgal/polygonal_surface_reconstruction.h>
//extras formal calculation
#include <Cgal/compute_average_spacing.h>
#include <Cgal/pca_estimate_normals.h>
#include <Cgal/mst_orient_normals.h>
#include <utility> // defines std::pair
#include <list>

///////////// TRY TO COMBINE norMAL CALculaTION + RANSAC + polyFIT

#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,int> PNI;
typedef std::vector<PNI> Point_vector;

//bij normals calc: Point with normal vector stored in a std::pair.
//typedef std::pair<Point,Vector> PointVectorPair;

typedef Cgal::Nth_of_tuple_property_map<0,PNI> Plane_index_map;

typedef std::pair<Point,Vector> PointVectorPair;

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::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;


// Concurrency
typedef Cgal::Parallel_if_available_tag Concurrency_tag;


std::vector<PNI> convert_function(const std::list<PointVectorPair>& list)
{
    std::vector<PNI> out;//a temporary object to store the output    
    for (std::size_t i = 0; i < list.size(); ++i) {
    {   
        out[i].get<0>() = list[i].get<0>();
        out[i].get<1>() = list[i].get<1>();
        
    }
    return out;
}

/*
* This example first extracts planes from the input point cloud
* (using RANSAC with default parameters) and then reconstructs
* the surface model from the planes.
*/

int main(int argc,char* argv[])
{
    const char* fname = (argc > 1) ? argv[1] : "data/47038.xyz";
    // Reads a .xyz point set file in points[].
    std::list<PointVectorPair> points;
    std::ifstream stream(fname);
    if (!stream ||
        !Cgal::read_xyz_points(stream,std::back_inserter(points),Cgal::parameters::point_map(Cgal::First_of_pair_property_map<PointVectorPair>())))
    {
        std::cerr << "Error: cannot read file " << fname << std::endl;
        return EXIT_FAILURE;
    }
    std::cout << "Calculating normals...";

    Cgal::Timer t;
    t.start();

    // Estimates normals direction.
    // Note: pca_estimate_normals() requires a range of points
    // as well as property maps to access each point's position and normal.
    const int nb_neighbors = 18; // K-nearest neighbors = 3 rings
    if (argc > 2 && std::strcmp(argv[2],"-r") == 0) // Use a fixed neighborhood radius
    {
        // First compute a spacing using the K parameter
        double spacing
            = Cgal::compute_average_spacing<Concurrency_tag>
            (points,nb_neighbors,Cgal::parameters::point_map(Cgal::First_of_pair_property_map<PointVectorPair>()));
        // Then,estimate normals with a fixed radius
        Cgal::pca_estimate_normals<Concurrency_tag>
            (points,// when using a neighborhood radius,K=0 means no limit on the number of neighbors returns
                Cgal::parameters::point_map(Cgal::First_of_pair_property_map<PointVectorPair>()).
                normal_map(Cgal::Second_of_pair_property_map<PointVectorPair>()).
                neighbor_radius(2. * spacing)); // use 2*spacing as neighborhood radius
    }
    else // Use a fixed number of neighbors
    {
        Cgal::pca_estimate_normals<Concurrency_tag>
            (points,Cgal::parameters::point_map(Cgal::First_of_pair_property_map<PointVectorPair>()).
                normal_map(Cgal::Second_of_pair_property_map<PointVectorPair>()));
    }

    // Orients normals.
    // Note: mst_orient_normals() requires a range of points
    // as well as property maps to access each point's position and normal.
    std::list<PointVectorPair>::iterator unoriented_points_begin =
        Cgal::mst_orient_normals(points,Cgal::parameters::point_map(Cgal::First_of_pair_property_map<PointVectorPair>()).
            normal_map(Cgal::Second_of_pair_property_map<PointVectorPair>()));
    // Optional: delete points with an unoriented normal
    // if you plan to call a reconstruction algorithm that expects oriented normals.
    points.erase(unoriented_points_begin,points.end());

    /// Convert pointvectorpair to pointvector

    
    std::vector<PNI> points2; // store points
    points2 = convert_function(points);
    ///
    /// use points + normals for surface reconstruction
    ///
    
    // Shape detection
    Efficient_ransac ransac;
    ransac.set_input(points2);
    ransac.add_shape_factory<Plane>();

    std::cout << "Extracting planes...";
    t.reset();
    ransac.detect();

    Efficient_ransac::Plane_range planes = ransac.planes();
    std::size_t num_planes = planes.size();

    std::cout << " Done. " << num_planes << " planes extracted. Time: " << t.time() << " sec." << std::endl;

    // Stores the plane index of each point as the third element of the tuple.
    Point_to_shape_index_map shape_index_map(points,planes);
    for (std::size_t i = 0; i < points.size(); ++i) {
        // Uses the get function from the property map that accesses the 3rd element of the tuple.
        int plane_index = get(shape_index_map,i);
        points2[i].get<2>() = plane_index;
    }

    //////////////////////////////////////////////////////////////////////////

    std::cout << "Generating candidate faces...";
    t.reset();

    polygonal_surface_reconstruction algo(
        points,Point_map(),normal_map(),Plane_index_map()
    );

    std::cout << " Done. Time: " << t.time() << " sec." << std::endl;

    //////////////////////////////////////////////////////////////////////////

    Surface_mesh model;

    std::cout << "Reconstructing...";
    t.reset();

    if (!algo.reconstruct<MIP_Solver>(model)) {
        std::cerr << " Failed: " << algo.error_message() << std::endl;
        return EXIT_FAILURE;
    }

    const std::string& output_file("data/cube_result_47038.off");
    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 EXIT_FAILURE;
    }

    //////////////////////////////////////////////////////////////////////////

    // Also stores the candidate faces as a surface mesh to a file
    Surface_mesh candidate_faces;
    algo.output_candidate_faces(candidate_faces);
    const std::string& candidate_faces_file("data/cube_candidate_faces_47038.off");
    std::ofstream candidate_stream(candidate_faces_file.c_str());
    if (candidate_stream && Cgal::write_off(candidate_stream,candidate_faces)) {
        // flush the buffer
        output_stream << std::flush;
        std::cout << "Candidate faces saved to " << candidate_faces_file << "." << std::endl;
    }

    return EXIT_SUCCESS;
}



#else

int main(int,char**)
{
    std::cerr << "This test requires either GLPK or SCIP.\n";
    return EXIT_SUCCESS;
}

#endif  // defined(Cgal_USE_GLPK) || defined(Cgal_USE_SCIP)

谢谢，

芭芭拉

解决方法

您尝试在列表中使用 operator[]：

std::vector<PNI> convert_function(const std::list<PointVectorPair>& list)
{
...
    for (std::size_t i = 0; i < list.size(); ++i)
    {   
        out[i].get<0>() = list[i].get<0>();

你根本不能做list[i]，因为std::list不提供它。

https://en.cppreference.com/w/cpp/container/list

相反，您可以使用：

auto it = list.begin();

和

it != list.end();

和

it++;

结合

(*it).get<0>();

允许您从头开始，检测结束并推进迭代器。结合对迭代器的取消引用。

在尝试写入之前，您还需要 resize() out 向量。