CGAL 3D Convex Hull с собственной реализацией «Черты» содержит внутренние точки - PullRequest
0 голосов
/ 23 октября 2018

Моя цель - вычислить с помощью CGAL статический выпуклый корпус из набора трехмерных точек.Я прочитал статью на правильной странице документации CGAL и реализовал первый алгоритм, в котором я получаю Polyhedron_3 (выпуклый корпус) из набора трехмерных точек, которые я ранее преобразовал в Kernel::Point_3 структуру.Это работает.

Затем я попытался применить алгоритм непосредственно к исходным точкам, , т.е. , без применения преобразования к структуре Kernel::Point_3.Из вышеупомянутой документации я понял, что способ сделать это - реализовать новый класс Черт, который определяет геометрические примитивы, используемые алгоритмом.Поэтому я следовал ConvexHullTraits_3 Ссылка на концепцию и, взглянув на файл Convex_hull_traits_3.h CGAL, я создал следующий код:

#include <iostream>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Convex_hull_traits_3.h>
#include <CGAL/internal/Projection_traits_3.h>
#include <CGAL/predicates/kernel_ftC3.h>
#include <CGAL/Kernel/global_functions_2.h>
#include <CGAL/Surface_mesh.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel My_kernel;

class My_point_3
{
private:

    typedef My_kernel::Point_3 Point_3;
    Point_3 point; //!< coordinates of the point
    int index; //!< additional info about the point

public:

    //constructors
    My_point_3() : index(-1), point() {};
    My_point_3(const double x, const double y, const double z, int i = -1) : index(i), point(x, y, z) {};
    My_point_3(Point_3& point, int i = -1) : index(i), point(point) {};
    My_point_3(Point_3 point, int i = -1) : index(i), point(point) {};

    //access to members for reading
    const int & Index() const { return index; }
    const Point_3 & Point()const { return point; }

    //access to members for writing
    int & Index() { return index; }
    Point_3 & Point() { return point; }

    //operators
    bool operator==(const My_point_3 &p) const
    {
        return (point == p.Point());
    }
    bool operator!=(const My_point_3 &p) const
    {
        return !(*this == p);
    }

    //see §5.3 of https://doc.cgal.org/latest/Kernel_23/index.html#sectionextensiblekernel
    const double & x() const { return point.x(); }
    const double & y() const { return point.y(); }
    const double & z() const { return point.z(); }
};

//projection traits implemented to satisfy https://doc.cgal.org/latest/Convex_hull_2/group__PkgConvexHull2Functions.html#gac13b4efbc337c7a8d5ad418521edcd4f
template <int dim>
class My_projection_traits_3
{
public:

    //Point_2
    typedef My_point_3                                Point_2;

    //Less_signed_distance_to_line_2
    typedef CGAL::internal::Less_signed_distance_to_line_projected_3<My_kernel, dim> Default_less_signed_distance_to_line_2;
    template <int d>
    class Less_signed_distance_to_line_projected_3
    {
    public:
        typedef My_point_3           Point_3;
        typedef My_kernel::Point_2   Point_2;
        My_kernel::FT x(const My_kernel::Point_3 &p) const { return CGAL::internal::Projector<My_kernel, d>::x(p); }
        My_kernel::FT y(const My_kernel::Point_3 &p) const { return CGAL::internal::Projector<My_kernel, d>::y(p); }

        Point_2 project(const Point_3& p) const
        {
            return Point_2(x(p.Point()), y(p.Point()));
        }

        bool operator()(const Point_3& p, const Point_3& q, const Point_3& r, const Point_3& s) const
        {
            Default_less_signed_distance_to_line_2 predicate;
            return predicate(project(p), project(q), project(r), project(s));
        }

        typedef bool result_type;
    };
    typedef Less_signed_distance_to_line_projected_3<dim>    Less_signed_distance_to_line_2;

    Less_signed_distance_to_line_2 less_signed_distance_to_line_2_object() const
    {
        return Less_signed_distance_to_line_2();
    }

    //Orientation_2
    template <int d>
    class Orientation_projected_3
    {
    public:
        typedef My_point_3           Point_3;
        typedef My_kernel::Point_2   Point_2;
        My_kernel::FT x(const My_kernel::Point_3 &p) const { return CGAL::internal::Projector<My_kernel, d>::x(p); }
        My_kernel::FT y(const My_kernel::Point_3 &p) const { return CGAL::internal::Projector<My_kernel, d>::y(p); }

        Point_2 project(const Point_3& p) const
        {
            return Point_2(x(p.Point()), y(p.Point()));
        }

        CGAL::Orientation operator()(const Point_3& p, const Point_3& q, const Point_3& r) const
        {
            return CGAL::orientation(project(p), project(q), project(r));
        }

        typedef CGAL::Orientation result_type;
    };
    typedef Orientation_projected_3<dim>                     Orientation_2;

    Orientation_2 orientation_2_object() const
    {
        return Orientation_2();
    }

    //Left_turn_2
    struct Left_turn_2
    {
    public:

        bool operator()(const Point_2& p, const Point_2& q, const Point_2& r) const
        {
            Orientation_2 predicate;
            return (predicate(p, q, r) == CGAL::LEFT_TURN);
        }

        typedef bool result_type;
    };

    Left_turn_2 left_turn_2_object() const
    {
        return Left_turn_2();
    }

    //Less_xy_2
    typedef typename CGAL::internal::Projector<My_kernel, dim>::Compare_x_2              Compare_x_2;
    typedef typename CGAL::internal::Projector<My_kernel, dim>::Less_y_2                 Less_y_2;
    struct Less_xy_2
    {
    public:

        bool operator()(const Point_2& p, const Point_2& q) const
        {
            Compare_x_2 predicate1;
            CGAL::Comparison_result crx = predicate1(p.Point(), q.Point());
            if (crx == CGAL::SMALLER) { return true; }
            if (crx == CGAL::LARGER) { return false; }
            Less_y_2 predicate2;
            return predicate2(p.Point(), q.Point());
        }

        typedef bool result_type;
    };

    Less_xy_2 less_xy_2_object() const
    {
        return Less_xy_2();
    }

    //Less_yx_2
    typedef typename CGAL::internal::Projector<My_kernel, dim>::Compare_y_2              Compare_y_2;
    typedef typename CGAL::internal::Projector<My_kernel, dim>::Less_x_2                 Less_x_2;
    struct Less_yx_2
    {
    public:

        bool operator()(const Point_2& p, const Point_2& q) const
        {
            Compare_y_2 predicate1;
            CGAL::Comparison_result cry = predicate1(p.Point(), q.Point());
            if (cry == CGAL::SMALLER) { return true; }
            if (cry == CGAL::LARGER) { return false; }
            Less_x_2 predicate2;
            return predicate2(p.Point(), q.Point());
        }

        typedef bool result_type;
    };

    Less_yx_2 less_yx_2_object() const
    {
        return Less_yx_2();
    }

    //Equal_2
    typedef typename CGAL::internal::Projector<My_kernel, dim>::Equal_x_2                Equal_x_2;
    typedef typename CGAL::internal::Projector<My_kernel, dim>::Equal_y_2                Equal_y_2;
    struct Equal_2
    {
    public:

        bool operator()(const Point_2& p, const Point_2& q) const
        {
            Equal_x_2 predicate1;
            Equal_y_2 predicate2;
            return predicate1(p.Point(), q.Point()) & predicate2(p.Point(), q.Point());
        }

        typedef bool result_type;
    };

    Equal_2 equal_2_object() const
    {
        return Equal_2();
    }
};

class Point_triple
{
protected:
    typedef My_kernel::FT                   FT;
    typedef My_kernel::Point_3              Point_3;
    typedef My_kernel::Vector_3             Vector_3;

    Point_3 p_, q_, r_;
public:

    Point_triple() {}

    Point_triple(const Point_3 &p, const Point_3 &q, const Point_3 &r): p_(p), q_(q), r_(r) {}

    const Point_3& p() const { return p_; }
    const Point_3& q() const { return q_; }
    const Point_3& r() const { return r_; }
};

typedef CGAL::Convex_hull_traits_3<My_kernel> Default_traits;

class My_traits
{
public:
    typedef My_point_3                    Point_3;
    typedef My_kernel::Segment_3          Segment_3;
    typedef My_kernel::Triangle_3         Triangle_3;
    typedef Point_triple                  Plane_3;

    class Construct_segment_3
    {
    public:

        Segment_3 operator ()(const Point_3& p, const Point_3& q) const
        {
            Default_traits::Construct_segment_3 object;
            return object(p.Point(), q.Point());
        }

        typedef Segment_3 result_type;
    };

    class Construct_plane_3
    {
    public:

        Plane_3 operator ()(const Point_3& p, const Point_3& q, const Point_3& r) const
        {
            CGAL::Convex_hull_traits_3<My_kernel>::Construct_plane_3 object;
            CGAL::Point_triple<My_kernel> plane = object(p.Point(), q.Point(), r.Point());
            return Point_triple(plane.p(), plane.q(), plane.r());
        }

        typedef Plane_3 result_type;
    };

    class Construct_triangle_3
    {
    public:

        Triangle_3 operator ()(const Point_3& p, const Point_3& q, const Point_3& r) const
        {
            Default_traits::Construct_triangle_3 object;
            return object(p.Point(), q.Point(), r.Point());
        }

        typedef Triangle_3 result_type;
    };

    class Equal_3
    {
    public:

        bool operator ()(const Point_3& p, const Point_3& q) const
        {
            Default_traits::Equal_3 predicate;
            return predicate(p.Point(), q.Point());
        }

        typedef bool result_type;
    };

    class Collinear_3
    {
    public:

        bool operator ()(const Point_3& p, const Point_3& q, const Point_3& r) const
        {
            Default_traits::Collinear_3 predicate;
            return predicate(p.Point(), q.Point(), r.Point());
        }

        typedef bool result_type;
    };

    class Coplanar_3
    {
    public:

        bool operator ()(const Point_3& p, const Point_3& q, const Point_3& r, const Point_3& s) const
        {
            Default_traits::Coplanar_3 predicate;
            return predicate(p.Point(), q.Point(), r.Point(), s.Point());
        }

        typedef bool result_type;
    };

    class Has_on_positive_side_3
    {
    public:

        bool operator()(const Plane_3& pl, const Point_3& p) const
        {
            Default_traits::Orientation_3 predicate;
            return (predicate(pl.p(), pl.q(), pl.r(), p.Point()) == CGAL::POSITIVE);
        }

        typedef bool result_type;
    };

    class Less_distance_to_point_3
    {
    public:

        bool operator ()(const Point_3 & p, const Point_3 & q, const Point_3 & r) const
        {
            Default_traits::Less_distance_to_point_3 predicate;
            return predicate(p.Point(), q.Point(), r.Point());
        }

        typedef bool result_type;
    };

    class Less_signed_distance_to_plane_3
    {

    public:

        bool operator()(const Plane_3 & h, const Point_3 & p, const Point_3 & q) const
        {
            const My_kernel::Point_3 & hp = h.p();
            const My_kernel::Point_3 & hq = h.q();
            const My_kernel::Point_3 & hr = h.r();
            return CGAL::has_smaller_signed_dist_to_planeC3(
                hp.x(), hp.y(), hp.z(), 
                hq.x(), hq.y(), hq.z(), 
                hr.x(), hr.y(), hr.z(), 
                p.Point().x(), p.Point().y(), 
                p.Point().z(), q.Point().x(), 
                q.Point().y(), q.Point().z()
            );
        }

        typedef bool result_type;
    };

    typedef My_projection_traits_3<2> Traits_xy_3;
    typedef My_projection_traits_3<0> Traits_yz_3;
    typedef My_projection_traits_3<1> Traits_xz_3;

    Construct_segment_3 construct_segment_3_object() const
    {
        return Construct_segment_3();
    }

    Construct_plane_3 construct_plane_3_object() const
    {
        return Construct_plane_3();
    }

    Construct_triangle_3 construct_triangle_3_object() const
    {
        return Construct_triangle_3();
    }

    Collinear_3 collinear_3_object() const
    {
        return Collinear_3();
    }

    Coplanar_3 coplanar_3_object() const
    {
        return Coplanar_3();
    }

    Less_distance_to_point_3 less_distance_to_point_3_object() const
    {
        return Less_distance_to_point_3();
    }

    Has_on_positive_side_3 has_on_positive_side_3_object() const
    {
        return Has_on_positive_side_3();
    }

    Equal_3 equal_3_object() const
    {
        return Equal_3();
    }

    Less_signed_distance_to_plane_3 less_signed_distance_to_plane_3_object() const
    {
        return Less_signed_distance_to_plane_3();
    }

    Traits_xy_3 construct_traits_xy_3_object() const
    {
        return Traits_xy_3();
    }

    Traits_yz_3 construct_traits_yz_3_object() const
    {
        return Traits_yz_3();
    }

    Traits_xz_3 construct_traits_xz_3_object() const
    {
        return Traits_xz_3();
    }
};


typedef My_kernel::Point_3 Point_3;
typedef CGAL::Surface_mesh<My_point_3>               Surface_mesh;

int main()
{
    std::vector<My_point_3> input_points;
    input_points.push_back(My_point_3(Point_3(0, 0, 0), 1));
    input_points.push_back(My_point_3(Point_3(1, 0, 0), 2));
    input_points.push_back(My_point_3(Point_3(0, 1, 0), 3));
    input_points.push_back(My_point_3(Point_3(0, 0, 1), 4));
    input_points.push_back(My_point_3(Point_3(0.5, 0.5, 0.5), 5));

    Surface_mesh mesh;
    CGAL::convex_hull_3(input_points.begin(), input_points.end(), mesh, My_traits());

    Surface_mesh::Vertex_iterator v;
    for (v = mesh.vertices_begin(); v != mesh.vertices_end(); v++)
    {
        std::cout << input_points[*v].Index() << " (" << input_points[*v].Point() << ")" << std::endl;
    }

    return 0;
}

Однако при попытке компиляции выдается следующая ошибка: [...] ( Редактировать : исправлено путем обновления CGAL с 4.12 до 4.13и некоторые другие исправления - см. историю изменений и комментарии ниже для получения подробной информации) .

Код компилируется;однако, когда я выполняю итерации по вершинам выпуклой оболочки в минимальном примере (см. main function выше), также печатаются вершины, не являющиеся выпуклой оболочкой.

Вопрос: чтоя делаю не так с моим кодом?чего не хватает?

...