surface_mesh_simplification.h

/********************************************************************
 * Copyright (C) 2020-2021 by Liangliang Nan <liangliang.nan@gmail.com>
 * Copyright (C) 2011-2020 the Polygon Mesh Processing Library developers.
 *
 * The code in this file is adapted from the PMP (Polygon Mesh Processing
 * Library) with modifications.
 *      https://github.com/pmp-library/pmp-library
 * The original code was distributed under a MIT-style license, see
 *      https://github.com/pmp-library/pmp-library/blob/master/LICENSE.txt
 ********************************************************************/
 
#ifndef EASY3D_ALGO_SURFACE_MESH_SIMPLIFICATION_H
#define EASY3D_ALGO_SURFACE_MESH_SIMPLIFICATION_H
 
#include <vector>
 
#include <easy3d/core/surface_mesh.h>
#include <easy3d/core/heap.h>
 
 
 
 
namespace easy3d {

    class SurfaceMeshSimplification {
    public:
        explicit SurfaceMeshSimplification(SurfaceMesh *mesh);

        ~SurfaceMeshSimplification();

        void initialize(float aspect_ratio = 0.0, float edge_length = 0.0,
                        unsigned int max_valence = 0, float normal_deviation = 0.0,
                        float hausdorff_error = 0.0);

        void simplify(unsigned int n_vertices);
 
    private:
        /*
                    vl
                    *
                   / \
                  /   \
                 / fl  \
             v0 *------>* v1
                 \ fr  /
                  \   /
                   \ /
                    *
                    vr
        */
        struct CollapseData {
            CollapseData(SurfaceMesh *sm, SurfaceMesh::Halfedge h);
 
            SurfaceMesh *mesh;
 
            SurfaceMesh::Halfedge v0v1; // Halfedge to be collapsed
            SurfaceMesh::Halfedge v1v0; // Reverse halfedge
            SurfaceMesh::Vertex v0;     // Vertex to be removed
            SurfaceMesh::Vertex v1;     // Remaining vertex
            SurfaceMesh::Face fl;       // Left face
            SurfaceMesh::Face fr;       // Right face
            SurfaceMesh::Vertex vl;     // Left vertex
            SurfaceMesh::Vertex vr;     // Right vertex
            SurfaceMesh::Halfedge v1vl, vlv0, v0vr, vrv1;
        };

        class HeapInterface {
        public:
            HeapInterface(SurfaceMesh::VertexProperty<float> prio, SurfaceMesh::VertexProperty<int> pos)
                    : prio_(prio), pos_(pos) {
            }
 
            bool less(SurfaceMesh::Vertex v0, SurfaceMesh::Vertex v1) { return prio_[v0] < prio_[v1]; }
 
            bool greater(SurfaceMesh::Vertex v0, SurfaceMesh::Vertex v1) { return prio_[v0] > prio_[v1]; }
 
            int get_heap_position(SurfaceMesh::Vertex v) { return pos_[v]; }
 
            void set_heap_position(SurfaceMesh::Vertex v, int pos) { pos_[v] = pos; }
 
        private:
            SurfaceMesh::VertexProperty<float> prio_;
            SurfaceMesh::VertexProperty<int> pos_;
        };

        class Quadric {
        public:
            Quadric(double a, double b, double c, double d,
                    double e, double f, double g,
                    double h, double i,
                    double j)
                    : a_(a), b_(b), c_(c), d_(d),
                      e_(e), f_(f), g_(g),
                      h_(h), i_(i),
                      j_(j) {}

            explicit Quadric(double a = 0.0, double b = 0.0, double c = 0.0, double d = 0.0)
                    : a_(a * a), b_(a * b), c_(a * c), d_(a * d),
                      e_(b * b), f_(b * c), g_(b * d),
                      h_(c * c), i_(c * d),
                      j_(d * d) {}

            Quadric(const vec3 &n, const vec3 &p) {
                *this = Quadric(n[0], n[1], n[2], -dot(n, p));
            }

            void clear() { a_ = b_ = c_ = d_ = e_ = f_ = g_ = h_ = i_ = j_ = 0.0; }

            Quadric &operator+=(const Quadric &q) {
                a_ += q.a_;
                b_ += q.b_;
                c_ += q.c_;
                d_ += q.d_;
                e_ += q.e_;
                f_ += q.f_;
                g_ += q.g_;
                h_ += q.h_;
                i_ += q.i_;
                j_ += q.j_;
                return *this;
            }

            Quadric &operator*=(double s) {
                a_ *= s;
                b_ *= s;
                c_ *= s;
                d_ *= s;
                e_ *= s;
                f_ *= s;
                g_ *= s;
                h_ *= s;
                i_ *= s;
                j_ *= s;
                return *this;
            }

            double operator()(const vec3 &p) const {
                const double x(p[0]), y(p[1]), z(p[2]);
                return a_ * x * x + 2.0 * b_ * x * y + 2.0 * c_ * x * z + 2.0 * d_ * x
                       + e_ * y * y + 2.0 * f_ * y * z + 2.0 * g_ * y
                       + h_ * z * z + 2.0 * i_ * z
                       + j_;
            }
 
        private:
 
            double a_, b_, c_, d_,
                    e_, f_, g_,
                    h_, i_,
                    j_;
        };
 
        //=============================================================================

        class NormalCone {
        public:
            NormalCone() = default;

            explicit NormalCone(const vec3 &normal, float angle = 0.0)
                    : center_normal_(normal), angle_(angle) {
            }

            const vec3 &center_normal() const { return center_normal_; }

            float angle() const { return angle_; }

            NormalCone &merge(const vec3 &n) { return merge(NormalCone(n)); }

            NormalCone &merge(const NormalCone &nc) {
                const float dp = dot(center_normal_, nc.center_normal_);
 
                // axes point in same direction
                if (dp > 0.99999) {
                    angle_ = std::max(angle_, nc.angle_);
                }
 
                    // axes point in opposite directions
                else if (dp < -0.99999) {
                    angle_ = static_cast<float>(2 * M_PI);
                } else {
                    // new angle
                    float center_angle = std::acos(dp);
                    float min_angle = std::min(-angle_, center_angle - nc.angle_);
                    float max_angle = std::max(angle_, center_angle + nc.angle_);
                    angle_ = 0.5f * (max_angle - min_angle);
 
                    // axis by SLERP
                    float axis_angle = 0.5f * (min_angle + max_angle);
                    center_normal_ = ((center_normal_ * std::sin(center_angle - axis_angle) +
                                       nc.center_normal_ * std::sin(axis_angle)) /
                                      std::sin(center_angle));
                }
 
                return *this;
            }
 
        private:
            vec3 center_normal_;
            float angle_;
        };
 
        typedef Heap<SurfaceMesh::Vertex, HeapInterface> PriorityQueue;
 
        typedef std::vector<vec3> Points;
 
    private:
        // put the vertex v in the priority queue
        void enqueue_vertex(SurfaceMesh::Vertex v);
 
        // is collapsing the halfedge h allowed?
        bool is_collapse_legal(const CollapseData &cd);
 
        // what is the priority of collapsing the halfedge h
        float priority(const CollapseData &cd);
 
        // postprocess halfedge collapse
        void postprocess_collapse(const CollapseData &cd);
 
        // compute aspect ratio for face f
        float aspect_ratio(SurfaceMesh::Face f) const;
 
        // compute distance from p to triangle f
        float distance(SurfaceMesh::Face f, const vec3 &p) const;
 
    private:
        SurfaceMesh *mesh_;
 
        bool initialized_;
 
        SurfaceMesh::VertexProperty<float> vpriority_;
        SurfaceMesh::VertexProperty<SurfaceMesh::Halfedge> vtarget_;
        SurfaceMesh::VertexProperty<int> heap_pos_;
        SurfaceMesh::VertexProperty<Quadric> vquadric_;
        SurfaceMesh::FaceProperty<NormalCone> normal_cone_;
        SurfaceMesh::FaceProperty<Points> face_points_;
 
        SurfaceMesh::VertexProperty<vec3> vpoint_;
        SurfaceMesh::FaceProperty<vec3> fnormal_;
        SurfaceMesh::VertexProperty<bool> vselected_;
        SurfaceMesh::VertexProperty<bool> vfeature_;
        SurfaceMesh::EdgeProperty<bool> efeature_;
 
        PriorityQueue *queue_;
 
        bool has_selection_;
        bool has_features_;
        float normal_deviation_;
        float hausdorff_error_;
        float aspect_ratio_;
        float edge_length_;
        unsigned int max_valence_;
    };
 
} // namespace easy3d
 
#endif  // EASY3D_ALGO_SURFACE_MESH_SIMPLIFICATION_H