PolyhedronShape_fwd.hpp

Go to the documentation of this file.

#ifndef MEDUSA_BITS_DOMAINS_POLYHEDRONSHAPE_FWD_HPP_
#define MEDUSA_BITS_DOMAINS_POLYHEDRONSHAPE_FWD_HPP_
 
#include <medusa/Config.hpp>
#include <medusa/bits/domains/discretization_helpers.hpp>
#include <medusa/bits/domains/discretization_helpers_advanced.hpp>
#include <medusa/bits/domains/GeneralFill.hpp>
#include <medusa/bits/domains/DomainShape.hpp>
#include <medusa/bits/domains/PolytopeShape.hpp>
#include <medusa/bits/types/Vec.hpp>
 
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/Side_of_triangle_mesh.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/compute_normal.h>
 
#include <fstream>
 
namespace mm {
 
template <typename vec_t>
class PolyhedronShape : public DomainShape<vec_t> {
    static_assert(vec_t::dim == 3, "Only available in 3 dimensions.");
    using base_t = DomainShape<vec_t>;  
 
    /* CGAL types */
    typedef CGAL::Simple_cartesian<double> K;  
    typedef K::Point_3 Point;  
    typedef K::Vector_3 Vector;  
    typedef CGAL::Surface_mesh<Point> Mesh;  
    typedef Mesh::Vertex_index VI;  
    typedef Mesh::Face_index FI;  
    typedef CGAL::AABB_face_graph_triangle_primitive<Mesh> Primitive;
    typedef CGAL::AABB_tree<CGAL::AABB_traits<K, Primitive>> Tree;
    typedef CGAL::Side_of_triangle_mesh<Mesh, K> PointInside;
 
    Mesh surface;  
    Tree tree;  
    PointInside inside_tester;  
 
    void init() {
        tree.insert(surface.faces_begin(), surface.faces_end(), surface);
        tree.accelerate_distance_queries();
 
        auto vnormals = surface.add_property_map<VI, Vector>("v:normals", CGAL::NULL_VECTOR).first;
        auto fnormals = surface.add_property_map<FI, Vector>("f:normals", CGAL::NULL_VECTOR).first;
        CGAL::Polygon_mesh_processing::compute_normals(surface, vnormals, fnormals);
    }
 
  public:
    using typename base_t::scalar_t;
    using typename base_t::vector_t;
    using base_t::dim;
    using base_t::discretizeBoundaryWithDensity;
    using base_t::discretizeBoundaryWithStep;
 
 
    explicit PolyhedronShape(const Mesh& surface)
            : surface(surface), tree(), inside_tester(tree) { init(); }
 
    PolyhedronShape(const PolyhedronShape& other)
            : surface(other.surface), tree(), inside_tester(tree) { init(); }
 
    PolyhedronShape& operator=(const PolyhedronShape& other) {
        surface = other.surface;
        tree.clear();
        init();
        return *this;
    }
 
    static PolyhedronShape fromOFF(const std::string& filename) {
        std::ifstream input(filename);
        assert_msg(input, "Failed opening file %s: %s", filename, strerror(errno));
        Mesh surface;
        input >> surface;
        assert_msg(!input.fail(), "Invalid OFF file - either parsing failed or the data "
                                  "does not represent a two-manifold.");
        assert_msg(input, "Failed reading file %s: %s", filename, strerror(errno));
        assert_msg(!surface.is_empty(), "Surface is empty.");
        assert_msg(CGAL::is_triangle_mesh(surface), "Surface must be a triangle mesh.");
        assert_msg(CGAL::is_closed(surface), "Surface must be closed (without holes).");
        return PolyhedronShape(surface);
    }
 
    bool contains(const vec_t& point) const override {
        return inside_tester(Point(point[0], point[1], point[2])) != CGAL::ON_UNBOUNDED_SIDE;
    }
 
    std::pair<vec_t, vec_t> bbox() const override {
        auto b = tree.bbox();
        return {{b.xmin(), b.ymin(), b.zmin()}, {b.xmax(), b.ymax(), b.zmax()}};
    }
 
    DomainDiscretization<vec_t> discretizeBoundaryWithDensity(
            const std::function<scalar_t(vec_t)>& dr, int type) const override;
 
    std::ostream& print(std::ostream& os) const override {
        return os << "PolyhedronShape with " << surface.num_vertices() << " vertices"
                  << " and " << surface.num_faces() << " faces.";
    }
 
    PolyhedronShape* clone() const override { return new PolyhedronShape(*this); }
 
    static int rgb2type(uint8_t r, uint8_t g, uint8_t b) {
        return -((static_cast<int>(r) << 16) + (static_cast<int>(g) << 8) + b) - 1;
    }
 
  private:
    vec_t getPoint(VI idx) const {
        auto p = surface.point(idx);
        return {p.x(), p.y(), p.z()};
    }
};
 
// This must be hidden during a doxygen run due to a bug in processing const declarations with
// differently named template parameters.
// See https://github.com/doxygen/doxygen/issues/8178
#ifndef DOXYGEN_RUN
 
template <typename Scalar>
class PolytopeShape<Vec<Scalar, 3>> : public PolyhedronShape<Vec<Scalar, 3>> {
    using base_t = PolyhedronShape<Vec<Scalar, 3>>;
    using typename base_t::PolyhedronShape;
  public:
    PolytopeShape(const base_t& shape) : base_t(shape) {}
};
#endif
 
 
}  // namespace mm
 
#endif  // MEDUSA_BITS_DOMAINS_POLYHEDRONSHAPE_FWD_HPP_