poly_mesh.h

/********************************************************************
 * Copyright (C) 2015 Liangliang Nan <liangliang.nan@gmail.com>
 * https://3d.bk.tudelft.nl/liangliang/
 *
 * This file is part of Easy3D. If it is useful in your research/work,
 * I would be grateful if you show your appreciation by citing it:
 * ------------------------------------------------------------------
 *      Liangliang Nan.
 *      Easy3D: a lightweight, easy-to-use, and efficient C++ library
 *      for processing and rendering 3D data.
 *      Journal of Open Source Software, 6(64), 3255, 2021.
 * ------------------------------------------------------------------
 *
 * Easy3D is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License Version 3
 * as published by the Free Software Foundation.
 *
 * Easy3D is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 ********************************************************************/
 
#ifndef EASY3D_CORE_POLYHEDRAL_MESH_H
#define EASY3D_CORE_POLYHEDRAL_MESH_H
 
#include <easy3d/core/model.h>
 
#include <set>
 
#include <easy3d/core/types.h>
#include <easy3d/core/property.h>
 
 
namespace easy3d {

    class PolyMesh : public Model
    {
 
    public: //------------------------------------------------------ topology types

        class BaseHandle
        {
        public:
            explicit BaseHandle(int _idx=-1) : idx_(_idx) {}

            int idx() const { return idx_; }

            void reset() { idx_=-1; }

            bool is_valid() const { return idx_ != -1; }

            bool operator==(const BaseHandle& _rhs) const {
                return idx_ == _rhs.idx_;
            }

            bool operator!=(const BaseHandle& _rhs) const {
                return idx_ != _rhs.idx_;
            }

            bool operator<(const BaseHandle& _rhs) const {
                return idx_ < _rhs.idx_;
            }

            struct Hash {
                std::size_t operator()(const BaseHandle& h) const { return h.idx(); }
            };
 
        private:
            friend class PolyMesh;
            int idx_;
        };

        struct Vertex : BaseHandle {
            explicit Vertex(int _idx=-1) : BaseHandle(_idx) {}
            std::ostream& operator<<(std::ostream& os) const { return os << 'v' << idx(); }
        };

        struct Edge : BaseHandle {
            explicit Edge(int _idx=-1) : BaseHandle(_idx) {}
            std::ostream& operator<<(std::ostream& os) const { return os << 'e' << idx(); }
        };

        struct HalfFace : BaseHandle {
            explicit HalfFace(int _idx=-1) : BaseHandle(_idx) {}
            std::ostream& operator<<(std::ostream& os) const { return os << 'h' << idx(); }
        };

        struct Face : BaseHandle {
            explicit Face(int _idx=-1) : BaseHandle(_idx) {}
            std::ostream& operator<<(std::ostream& os) const { return os << 'f' << idx(); }
        };

        struct Cell : BaseHandle {
            explicit Cell(int _idx=-1) : BaseHandle(_idx) {}
            std::ostream& operator<<(std::ostream& os) const { return os << 'c' << idx(); }
        };
 
 
    public: //-------------------------------------------------- connectivity types

        struct VertexConnectivity {
            std::set<Vertex>     vertices_;  
            std::set<Edge>       edges_;     
            std::set<HalfFace>   halffaces_; 
            std::set<Cell>       cells_;     
            void read(std::istream& in);
            void write(std::ostream& out) const;
        };

        struct EdgeConnectivity {
            std::vector<Vertex>  vertices_;  
            std::set<HalfFace>   halffaces_; 
            std::set<Cell>       cells_;     
            void read(std::istream& in);
            void write(std::ostream& out) const;
        };

        struct HalfFaceConnectivity {
            std::vector<Vertex> vertices_; 
            std::set<Edge>      edges_;    
            Cell                cell_;     
            HalfFace            opposite_; 
            void read(std::istream& in);
            void write(std::ostream& out) const;
        };

        struct CellConnectivity {
            std::set<Vertex>        vertices_;    
            std::set<Edge>          edges_;       
            std::vector<HalfFace>   halffaces_;   
            void read(std::istream& in);
            void write(std::ostream& out) const;
        };
 
 
    public: //------------------------------------------------------ property types

        template <class T> class VertexProperty : public Property<T>
        {
        public:
            VertexProperty() = default;
            explicit VertexProperty(Property<T> p) : Property<T>(p) {}

            typename Property<T>::reference operator[](Vertex v) {
                return Property<T>::operator[](v.idx());
            }

            typename Property<T>::const_reference operator[](Vertex v) const {
                return Property<T>::operator[](v.idx());
            }
        };

        template <class T> class EdgeProperty : public Property<T>
        {
        public:
            EdgeProperty() = default;
            explicit EdgeProperty(Property<T> p) : Property<T>(p) {}

            typename Property<T>::reference operator[](Edge e) {
                return Property<T>::operator[](e.idx());
            }

            typename Property<T>::const_reference operator[](Edge e) const {
                return Property<T>::operator[](e.idx());
            }
        };

        template <class T> class HalfFaceProperty : public Property<T>
        {
        public:
            HalfFaceProperty() = default;
            explicit HalfFaceProperty(Property<T> p) : Property<T>(p) {}

            typename Property<T>::reference operator[](HalfFace h) {
                return Property<T>::operator[](h.idx());
            }

            typename Property<T>::const_reference operator[](HalfFace h) const {
                return Property<T>::operator[](h.idx());
            }
        };

        template <class T> class FaceProperty : public Property<T>
        {
        public:
            FaceProperty() = default;
            explicit FaceProperty(Property<T> p) : Property<T>(p) {}

            typename Property<T>::reference operator[](Face f) {
                return Property<T>::operator[](f.idx());
            }

            typename Property<T>::const_reference operator[](Face f) const {
                return Property<T>::operator[](f.idx());
            }
        };

        template <class T> class CellProperty : public Property<T>
        {
        public:
            CellProperty() = default;
            explicit CellProperty(Property<T> p) : Property<T>(p) {}

            typename Property<T>::reference operator[](Cell c) {
                return Property<T>::operator[](c.idx());
            }

            typename Property<T>::const_reference operator[](Cell c) const {
                return Property<T>::operator[](c.idx());
            }
        };
 

        template <class T> class ModelProperty : public Property<T>
        {
        public:
            ModelProperty() = default;
            explicit ModelProperty(Property<T> p) : Property<T>(p) {}

            typename Property<T>::reference operator[](size_t idx) {
                return Property<T>::operator[](idx);
            }

            typename Property<T>::const_reference operator[](size_t idx) const {
                return Property<T>::operator[](idx);
            }
        };
 
 
    public: //------------------------------------------------------ iterator types
        class VertexIterator
        {
        public:
            explicit VertexIterator(Vertex v=Vertex(), const PolyMesh* m=nullptr) : hnd_(v), mesh_(m) { }

            Vertex operator*()  const { return  hnd_; }

            bool operator==(const VertexIterator& rhs) const {
                return (hnd_==rhs.hnd_);
            }

            bool operator!=(const VertexIterator& rhs) const {
                return !operator==(rhs);
            }

            VertexIterator& operator++() {
                ++hnd_.idx_;
                assert(mesh_);
                return *this;
            }

            VertexIterator& operator--() {
                --hnd_.idx_;
                assert(mesh_);
                return *this;
            }
 
        private:
            Vertex  hnd_;
            const PolyMesh* mesh_;
        };

        class EdgeIterator
        {
        public:
            explicit EdgeIterator(Edge e=Edge(), const PolyMesh* m=nullptr) : hnd_(e), mesh_(m) { }

            Edge operator*()  const { return  hnd_; }

            bool operator==(const EdgeIterator& rhs) const {
                return (hnd_==rhs.hnd_);
            }

            bool operator!=(const EdgeIterator& rhs) const {
                return !operator==(rhs);
            }

            EdgeIterator& operator++() {
                ++hnd_.idx_;
                assert(mesh_);
                return *this;
            }

            EdgeIterator& operator--() {
                --hnd_.idx_;
                assert(mesh_);
                return *this;
            }
 
        private:
            Edge  hnd_;
            const PolyMesh* mesh_;
        };
 

        class HalfFaceIterator
        {
        public:
            explicit HalfFaceIterator(HalfFace h=HalfFace(), const PolyMesh* m=nullptr) : hnd_(h), mesh_(m) { }

            HalfFace operator*()  const { return  hnd_; }

            bool operator==(const HalfFaceIterator& rhs) const {
                return (hnd_==rhs.hnd_);
            }

            bool operator!=(const HalfFaceIterator& rhs) const {
                return !operator==(rhs);
            }

            HalfFaceIterator& operator++() {
                ++hnd_.idx_;
                assert(mesh_);
                return *this;
            }

            HalfFaceIterator& operator--() {
                --hnd_.idx_;
                assert(mesh_);
                return *this;
            }
 
        private:
            HalfFace  hnd_;
            const PolyMesh* mesh_;
        };
 

        class FaceIterator
        {
        public:
            explicit FaceIterator(Face f=Face(), const PolyMesh* m=nullptr) : hnd_(f), mesh_(m) { }

            Face operator*()  const { return  hnd_; }

            bool operator==(const FaceIterator& rhs) const {
                return (hnd_==rhs.hnd_);
            }

            bool operator!=(const FaceIterator& rhs) const {
                return !operator==(rhs);
            }

            FaceIterator& operator++() {
                ++hnd_.idx_;
                assert(mesh_);
                return *this;
            }

            FaceIterator& operator--() {
                --hnd_.idx_;
                assert(mesh_);
                return *this;
            }
 
        private:
            Face  hnd_;
            const PolyMesh* mesh_;
        };
 

        class CellIterator
        {
        public:
            explicit CellIterator(Cell c=Cell(), const PolyMesh* m=nullptr) : hnd_(c), mesh_(m) { }

            Cell operator*() const { return  hnd_; }

            bool operator==(const CellIterator& rhs) const {
                return (hnd_==rhs.hnd_);
            }

            bool operator!=(const CellIterator& rhs) const {
                return !operator==(rhs);
            }

            CellIterator& operator++() {
                ++hnd_.idx_;
                assert(mesh_);
                return *this;
            }

            CellIterator& operator--() {
                --hnd_.idx_;
                assert(mesh_);
                return *this;
            }
 
        private:
            Cell  hnd_;
            const PolyMesh* mesh_;
        };
 
 
    public: //-------------------------- containers for C++11 range-based for loops

        class VertexContainer
        {
        public:
            VertexContainer(VertexIterator _begin, VertexIterator _end) : begin_(_begin), end_(_end) {}
            VertexIterator begin() const { return begin_; }
            VertexIterator end()   const { return end_;   }
        private:
            VertexIterator begin_, end_;
        };
 

        class EdgeContainer
        {
        public:
            EdgeContainer(EdgeIterator _begin, EdgeIterator _end) : begin_(_begin), end_(_end) {}
            EdgeIterator begin() const { return begin_; }
            EdgeIterator end()   const { return end_;   }
        private:
            EdgeIterator begin_, end_;
        };
 

        class HalffaceContainer
        {
        public:
            HalffaceContainer(HalfFaceIterator _begin, HalfFaceIterator _end) : begin_(_begin), end_(_end) {}
            HalfFaceIterator begin() const { return begin_; }
            HalfFaceIterator end()   const { return end_;   }
        private:
            HalfFaceIterator begin_, end_;
        };
 

        class FaceContainer
        {
        public:
            FaceContainer(FaceIterator _begin, FaceIterator _end) : begin_(_begin), end_(_end) {}
            FaceIterator begin() const { return begin_; }
            FaceIterator end()   const { return end_;   }
        private:
            FaceIterator begin_, end_;
        };
 
        
        class CellContainer
        {
        public:
            CellContainer(CellIterator _begin, CellIterator _end) : begin_(_begin), end_(_end) {}
            CellIterator begin() const { return begin_; }
            CellIterator end()   const { return end_;   }
        private:
            CellIterator begin_, end_;
        };
 
 
    public: //-------------------------------------------- constructor / destructor



        PolyMesh();

        ~PolyMesh() override = default;

        PolyMesh(const PolyMesh& rhs) { operator=(rhs); }

        PolyMesh& operator=(const PolyMesh& rhs);

        PolyMesh& assign(const PolyMesh& rhs);


        bool read(const std::string& filename);
        bool write(const std::string& filename) const;
 
    public: //----------------------------------------------- add new vertex / face / cell



        Vertex add_vertex(const vec3& p) { Vertex v = new_vertex(); vpoint_[v] = p; return v; }
        Cell add_cell(const std::vector<HalfFace>& faces);
        Cell add_tetra(HalfFace f0, HalfFace f1, HalfFace f2, HalfFace f3) { return add_cell({f0, f1, f2, f3}); }
        Cell add_tetra(Vertex v0, Vertex v1, Vertex v2, Vertex v3);
        Cell add_hexa(Vertex v0, Vertex v1, Vertex v2, Vertex v3, Vertex v4, Vertex v5, Vertex v6, Vertex v7);
        HalfFace add_face(const std::vector<Vertex>& vertices);
        HalfFace add_triangle(Vertex v0, Vertex v1, Vertex v2) { return add_face({ v0, v1, v2 }); }
        HalfFace add_quad(Vertex v0, Vertex v1, Vertex v2, Vertex v3) { return add_face({ v0, v1, v2, v3}); }

 
    public: //--------------------------------------------------- memory management



        unsigned int n_vertices() const { return (unsigned int) vprops_.size(); }
        unsigned int n_edges() const { return (unsigned int) eprops_.size(); }
        unsigned int n_halffaces() const { return (unsigned int) hprops_.size(); }
        unsigned int n_faces() const { return (unsigned int) fprops_.size(); }
        unsigned int n_cells() const { return (unsigned int) cprops_.size(); }

        void clear();

        void resize(unsigned int nv, unsigned int ne, unsigned int nf, unsigned int nc) {
            vprops_.resize(nv);
            eprops_.resize(ne);
            hprops_.resize(2 * nf);
            fprops_.resize(nf);
            cprops_.resize(nc);
        }

        bool is_valid(Vertex v) const {
            return (0 <= v.idx()) && (v.idx() < (int)n_vertices());
        }

        bool is_valid(Edge e) const {
            return (0 <= e.idx()) && (e.idx() < (int)n_edges());
        }

        bool is_valid(HalfFace h) const {
            return (0 <= h.idx()) && (h.idx() < (int)n_halffaces());
        }

        bool is_valid(Face f) const {
            return (0 <= f.idx()) && (f.idx() < (int)n_faces());
        }

        bool is_valid(Cell c) const {
            return (0 <= c.idx()) && (c.idx() < (int)n_cells());
        }

        
 
    public: //--------------------------------------------------- property handling



        template <class T> VertexProperty<T> add_vertex_property(const std::string& name, const T t=T()) {
            return VertexProperty<T>(vprops_.add<T>(name, t));
        }

        template <class T> EdgeProperty<T> add_edge_property(const std::string& name, const T t=T()) {
            return EdgeProperty<T>(eprops_.add<T>(name, t));
        }

        template <class T> HalfFaceProperty<T> add_halfface_property(const std::string& name, const T t=T()) {
            return HalfFaceProperty<T>(hprops_.add<T>(name, t));
        }

        template <class T> FaceProperty<T> add_face_property(const std::string& name, const T t=T()) {
            return FaceProperty<T>(fprops_.add<T>(name, t));
        }

        template <class T> CellProperty<T> add_cell_property(const std::string& name, const T t=T()) {
            return CellProperty<T>(cprops_.add<T>(name, t));
        }

        template <class T> ModelProperty<T> add_model_property(const std::string& name, const T t=T()) {
            return ModelProperty<T>(mprops_.add<T>(name, t));
        }

        template <class T> VertexProperty<T> get_vertex_property(const std::string& name) const {
            return VertexProperty<T>(vprops_.get<T>(name));
        }

        template <class T> EdgeProperty<T> get_edge_property(const std::string& name) const {
            return EdgeProperty<T>(eprops_.get<T>(name));
        }

        template <class T> HalfFaceProperty<T> get_halfface_property(const std::string& name) const {
            return HalfFaceProperty<T>(hprops_.get<T>(name));
        }

        template <class T> FaceProperty<T> get_face_property(const std::string& name) const {
            return FaceProperty<T>(fprops_.get<T>(name));
        }

        template <class T> CellProperty<T> get_cell_property(const std::string& name) const {
            return CellProperty<T>(cprops_.get<T>(name));
        }

        template <class T> ModelProperty<T> get_model_property(const std::string& name) const {
            return ModelProperty<T>(mprops_.get<T>(name));
        }

        template <class T> VertexProperty<T> vertex_property(const std::string& name, const T t=T()) {
            return VertexProperty<T>(vprops_.get_or_add<T>(name, t));
        }

        template <class T> EdgeProperty<T> edge_property(const std::string& name, const T t=T()) {
            return EdgeProperty<T>(eprops_.get_or_add<T>(name, t));
        }

        template <class T> HalfFaceProperty<T> halfface_property(const std::string& name, const T t=T()) {
            return HalfFaceProperty<T>(hprops_.get_or_add<T>(name, t));
        }

        template <class T> FaceProperty<T> face_property(const std::string& name, const T t=T()) {
            return FaceProperty<T>(fprops_.get_or_add<T>(name, t));
        }

        template <class T> CellProperty<T> cell_property(const std::string& name, const T t=T()) {
            return CellProperty<T>(cprops_.get_or_add<T>(name, t));
        }

        template <class T> ModelProperty<T> model_property(const std::string& name, const T t=T()) {
            return ModelProperty<T>(mprops_.get_or_add<T>(name, t));
        }
 

        template<class T>
        bool remove_vertex_property(VertexProperty<T> &p) { return vprops_.remove(p); }
        bool remove_vertex_property(const std::string &n) { return vprops_.remove(n); }
        template<class T>
        bool remove_edge_property(EdgeProperty<T> &p) { return eprops_.remove(p); }
        bool remove_edge_property(const std::string &n) { return eprops_.remove(n); }
        template<class T>
        bool remove_halfface_property(HalfFaceProperty<T> &p) { return hprops_.remove(p); }
        bool remove_halfface_property(const std::string &n) { return hprops_.remove(n); }
        template<class T>
        bool remove_face_property(FaceProperty<T> &p) { return fprops_.remove(p); }
        bool remove_face_property(const std::string &n) { return fprops_.remove(n); }
        template<class T>
        bool remove_cell_property(CellProperty<T> &p) { return cprops_.remove(p); }
        bool remove_cell_property(const std::string &n) { return cprops_.remove(n); }
        template<class T>
        bool remove_model_property(ModelProperty<T> &p) { return mprops_.remove(p); }
        bool remove_model_property(const std::string &n) { return mprops_.remove(n); }

        bool rename_vertex_property(const std::string &old_name, const std::string &new_name) {
            return vprops_.rename(old_name, new_name);
        }

        bool rename_halfface_property(const std::string &old_name, const std::string &new_name) {
            return hprops_.rename(old_name, new_name);
        }

        bool rename_face_property(const std::string &old_name, const std::string &new_name) {
            return fprops_.rename(old_name, new_name);
        }

        bool rename_edge_property(const std::string &old_name, const std::string &new_name) {
            return eprops_.rename(old_name, new_name);
        }

        bool rename_cell_property(const std::string &old_name, const std::string &new_name) {
            return cprops_.rename(old_name, new_name);
        }

        bool rename_model_property(const std::string &old_name, const std::string &new_name) {
            return mprops_.rename(old_name, new_name);
        }

        const std::type_info& get_vertex_property_type(const std::string& name) const {
            return vprops_.get_type(name);
        }

        const std::type_info& get_edge_property_type(const std::string& name) const {
            return eprops_.get_type(name);
        }

        const std::type_info& get_halfface_property_type(const std::string& name) const {
            return hprops_.get_type(name);
        }

        const std::type_info& get_face_property_type(const std::string& name) const {
            return fprops_.get_type(name);
        }

        const std::type_info& get_cell_property_type(const std::string& name) const {
            return cprops_.get_type(name);
        }

        const std::type_info& get_model_property_type(const std::string& name) const {
            return mprops_.get_type(name);
        }
 

        std::vector<std::string> vertex_properties() const {
            return vprops_.properties();
        }

        std::vector<std::string> edge_properties() const {
            return eprops_.properties();
        }

        std::vector<std::string> halfface_properties() const {
            return hprops_.properties();
        }

        std::vector<std::string> face_properties() const {
            return fprops_.properties();
        }

        std::vector<std::string> cell_properties() const {
            return cprops_.properties();
        }

        std::vector<std::string> model_properties() const {
            return mprops_.properties();
        }

        void property_stats(std::ostream &output) const override;

 
 
    public: //--------------------------------------------- iterators & circulators



        VertexIterator vertices_begin() const {
            return VertexIterator(Vertex(0), this);
        }

        VertexIterator vertices_end() const {
            return VertexIterator(Vertex(static_cast<int>(n_vertices())), this);
        }

        VertexContainer vertices() const {
            return VertexContainer(vertices_begin(), vertices_end());
        }

        EdgeIterator edges_begin() const {
            return EdgeIterator(Edge(0), this);
        }

        EdgeIterator edges_end() const {
            return EdgeIterator(Edge(static_cast<int>(n_edges())), this);
        }

        EdgeContainer edges() const {
            return EdgeContainer(edges_begin(), edges_end());
        }

        HalfFaceIterator halffaces_begin() const {
            return HalfFaceIterator(HalfFace(0), this);
        }

        HalfFaceIterator halffaces_end() const {
            return HalfFaceIterator(HalfFace(static_cast<int>(n_halffaces())), this);
        }

        HalffaceContainer halffaces() const {
            return HalffaceContainer(halffaces_begin(), halffaces_end());
        }

        FaceIterator faces_begin() const {
            return FaceIterator(Face(0), this);
        }

        FaceIterator faces_end() const {
            return FaceIterator(Face(static_cast<int>(n_faces())), this);
        }

        FaceContainer faces() const {
            return FaceContainer(faces_begin(), faces_end());
        }

        CellIterator cells_begin() const {
            return CellIterator(Cell(0), this);
        }

        CellIterator cells_end() const {
            return CellIterator(Cell(static_cast<int>(n_cells())), this);
        }

        CellContainer cells() const {
            return CellContainer(cells_begin(), cells_end());
        }

 
    public: //--------------------------------------------- adjacency access



        const std::set<Vertex>& vertices(Vertex v) const {
            return vconn_[v].vertices_;
        }

        HalfFace halfface(Face f, unsigned int i) const {
            assert(i<=1);
            return HalfFace(static_cast<int>((f.idx() << 1) + i));
        }

        Face face(HalfFace h) const {
            return Face((h.idx() >> 1));
        }

        HalfFace opposite(HalfFace h) const {
            return hconn_[h].opposite_;
        }

        Vertex vertex(Edge e, unsigned int i) const {
            assert(i<=1);
            return econn_[e].vertices_[i];
        }

        const std::vector<Vertex>& vertices(HalfFace h) const {
            return hconn_[h].vertices_;
        }

        const std::vector<Vertex>& vertices(Face f) const {
            return vertices(halfface(f, 0));
        }

        const std::set<Vertex>& vertices(Cell c) const {
            return cconn_[c].vertices_;
        }

        const std::set<Edge>& edges(Vertex v) const {
            return vconn_[v].edges_;
        }

        const std::set<Edge>& edges(HalfFace h) const {
            return hconn_[h].edges_;
        }

        const std::set<Edge>& edges(Cell c) const {
            return cconn_[c].edges_;
        }
        
        const std::set<HalfFace>& halffaces(Vertex v) const {
            return vconn_[v].halffaces_;
        }

        const std::set<HalfFace>& halffaces(Edge e) const {
            return econn_[e].halffaces_;
        }

        const std::vector<HalfFace>& halffaces(Cell c) const {
            return cconn_[c].halffaces_;
        }

        const std::set<Cell>& cells(Vertex v) const {
            return vconn_[v].cells_;
        }

        const std::set<Cell>& cells(Edge e) const {
            return econn_[e].cells_;
        }

        Cell cell(HalfFace h) const {
            return hconn_[h].cell_;
        }

 
    public: //--------------------------------------------- higher-level operations



        bool is_tetraheral_mesh() const;

        bool is_border(Vertex v) const {
            for (auto h : halffaces(v)) { if (is_border(h)) return true; }
            return false;
        }

        bool is_border(Edge e) const {
            for (auto h : halffaces(e)) { if (is_border(h)) return true; }
            return false;
        }

        bool is_border(HalfFace h) const {
            return (!cell(h).is_valid());
        }

        bool is_border(Face f) const {
            return is_border(halfface(f, 0)) || is_border(halfface(f, 1));
        }

        Edge find_edge(Vertex a, Vertex b) const;

        HalfFace find_half_face(const std::vector<Vertex>& vertices) const;

        bool is_degenerate(Face f) const;

        void extract_boundary(std::vector< std::vector<Vertex> >& faces) const;
 
    public: //------------------------------------------ geometry-related functions



        const vec3& position(Vertex v) const { return vpoint_[v]; }

        const std::vector<vec3>& points() const override { return vpoint_.vector(); }

        std::vector<vec3>& points() override { return vpoint_.vector(); }

        void update_face_normals();

        vec3 compute_face_normal(Face f) const;

        void update_vertex_normals();

        float edge_length(Edge e) const;

 
 
    private: //---------------------------------------------- allocate new elements

        Vertex new_vertex()
        {
            vprops_.push_back();
            return Vertex(static_cast<int>(n_vertices()-1));
        }

        Edge new_edge(Vertex s, Vertex t)
        {
            assert(s != t);
            eprops_.push_back();
            Edge e = Edge(static_cast<int>(n_edges() - 1));
            econn_[e].vertices_ = {s, t};
            vconn_[s].edges_.insert(e);
            vconn_[t].edges_.insert(e);
            vconn_[s].vertices_.insert(t);
            vconn_[t].vertices_.insert(s);
            return e;
        }

        HalfFace new_face()
        {
            fprops_.push_back();
            hprops_.push_back();
            hprops_.push_back();
            HalfFace h0(static_cast<int>(n_halffaces()-2));
            HalfFace h1(static_cast<int>(n_halffaces()-1));
 
            hconn_[h0].opposite_ = h1;
            hconn_[h1].opposite_ = h0;
 
            return h0;
        }

        Cell new_cell()
        {
            cprops_.push_back();
            return Cell(static_cast<int>(n_cells()-1));
        }
 
    private: //------------------------------------------------------- private data
 
        PropertyContainer vprops_;
        PropertyContainer eprops_;
        PropertyContainer hprops_;
        PropertyContainer fprops_;
        PropertyContainer cprops_;
        PropertyContainer mprops_;
 
        VertexProperty<VertexConnectivity>      vconn_;
        EdgeProperty<EdgeConnectivity>          econn_;
        HalfFaceProperty<HalfFaceConnectivity>  hconn_;
        CellProperty<CellConnectivity>          cconn_;
 
        VertexProperty<vec3>    vpoint_;
    };
 
 
    //------------------------------------------------------------ output operators

    inline std::ostream& operator<<(std::ostream& os, PolyMesh::Vertex v)
    {
        return (os << 'v' << v.idx());
    }

    inline std::ostream& operator<<(std::ostream& os, PolyMesh::Edge e)
    {
        return (os << 'e' << e.idx());
    }

    inline std::ostream& operator<<(std::ostream& os, PolyMesh::HalfFace h)
    {
        return (os << 'h' << h.idx());
    }

    inline std::ostream& operator<<(std::ostream& os, PolyMesh::Face f)
    {
        return (os << 'f' << f.idx());
    }

    inline std::ostream& operator<<(std::ostream& os, PolyMesh::Cell c)
    {
        return (os << 'c' << c.idx());
    }
 
} // namespace easy3d
 
#endif // EASY3D_CORE_POLYHEDRAL_MESH_H