首页 > 解决方案 > 给定点,计算法线并在 CGAL 中执行多边形表面重建

问题描述

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

但是,两者的输出和输入类型不同:

正常估计的输出是向量对的列表,参见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, Vector, 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> 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 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,
                0, // 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, nb_neighbors,
                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, nb_neighbors,
            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)

谢谢,

芭芭拉

标签: c++type-conversioncgalnormals

解决方案


您尝试在列表中使用 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


推荐阅读