// ========================================================================== //
// Class inclusion
#include "Utilities/mesh_operations.hpp"

// Local libraries
#include "Utilities/math_operations.hpp"

// External libraries
#include <GEL/CGLA/CGLA.h>
#include <GEL/HMesh/HMesh.h>

// Standard libraries
#include <queue>
#include <vector>

// Parallelisation
#include <omp.h>
#if USE_MPI == 1
    #include <mpi.h>
#endif

// Macros
#ifndef MPI_VERSION
    #define MPI_VERSION 0
#endif

namespace Utilities {

    // ====================================================================== //
    // Additional operations to work on meshes.
    int remove_small_components(HMesh::Manifold &m, int comp_size_thresh) {
        int work = 0;

        // Assign a number to all faces such that two faces are given the same
        // number if it is possible to reach one from the other via a path that
        // consists of faces which are connected by shared edges. For each
        // segment (group of faces with the same number) we count the number of
        // faces that belong to it.
        HMesh::FaceAttributeVector<int> face_segment(m.allocated_faces(), -1);
        int                             seg_no = 0;
        std::vector<int>                segment_sizes(0);
        for (auto fid : m.faces()) {
            if (face_segment[fid] == -1) {
                std::queue<HMesh::FaceID> fq;
                fq.push(fid);
                segment_sizes.push_back(1);
                face_segment[fid] = seg_no;
                while (!fq.empty()) {
                    HMesh::FaceID f = fq.front();
                    fq.pop();
                    circulate_face_ccw(m, f, [&](HMesh::FaceID fn) {
                        if (fn != HMesh::InvalidFaceID &&
                            face_segment[fn] == -1) {
                            face_segment[fn] = seg_no;
                            fq.push(fn);
                            ++segment_sizes[seg_no];
                        }
                    });
                }
                ++seg_no;
            }
        }

        // Remove all faces belonging to segments with a number of faces smaller
        // than the threshold.
        bool unclean = false;
        for (auto f : m.faces())
            if (segment_sizes[face_segment[f]] < comp_size_thresh) {
                m.remove_face(f);
                ++work;
                unclean = true;
            }
        if (unclean)
            m.cleanup(); // remove deleted entities from data structures
        return work;
    }

    int keep_largest_component(HMesh::Manifold &m) {
        int work = 0;

        // Assign a number to all faces such that two faces are given the same
        // number if it is possible to reach one from the other via a path that
        // consists of faces which are connected by shared edges. For each
        // segment (group of faces with the same number) we count the number of
        // faces that belong to it.
        HMesh::FaceAttributeVector<int> face_segment(m.allocated_faces(), -1);
        int                             seg_no = 0;
        std::vector<int>                segment_sizes(0);
        for (auto fid : m.faces()) {
            if (face_segment[fid] == -1) {
                std::queue<HMesh::FaceID> fq;
                fq.push(fid);
                segment_sizes.push_back(1);
                face_segment[fid] = seg_no;
                while (!fq.empty()) {
                    HMesh::FaceID f = fq.front();
                    fq.pop();
                    circulate_face_ccw(m, f, [&](HMesh::FaceID fn) {
                        if (fn != HMesh::InvalidFaceID &&
                            face_segment[fn] == -1) {
                            face_segment[fn] = seg_no;
                            fq.push(fn);
                            ++segment_sizes[seg_no];
                        }
                    });
                }
                ++seg_no;
            }
        }

        int comp_size_thresh = max(segment_sizes);

        // Remove all faces belonging to segments with a number of faces smaller
        // than the threshold.
        bool unclean = false;
        for (auto f : m.faces())
            if (segment_sizes[face_segment[f]] < comp_size_thresh) {
                m.remove_face(f);
                ++work;
                unclean = true;
            }
        if (unclean)
            m.cleanup(); // remove deleted entities from data structures
        return work;
    }

    void smooth_boundary(HMesh::Manifold &m, double angle_threshold,
                         int maxiter) {

        if (m.no_faces() <= 1) return;

        double avg = average_edge_length(m);

        int  iter   = 0;
        bool change = true;
        while (change && iter++ < maxiter) {
            change = false;

            // Initialize the pairing attribute std::vector
            int                               pair_value = 0;
            HMesh::VertexAttributeVector<int> pair(m.allocated_vertices());
            HMesh::VertexAttributeVector<int> visited(m.allocated_vertices(),
                                                      0);
            for (HMesh::VertexID vid : m.vertices()) pair[vid] = pair_value++;

            // Find non-convex regions on the boundary
            for (HMesh::HalfEdgeID hid : m.halfedges()) {
                if (!boundary(m, hid)) continue;

                // Make sure we have the outside halfedge
                HMesh::Walker w = m.walker(hid);
                if (w.face() != HMesh::InvalidFaceID) w = w.opp();

                // Fetch vertices, points and directions
                HMesh::VertexID vid1 = w.opp().vertex();
                HMesh::VertexID vid2 = w.vertex();
                HMesh::VertexID vid3 = w.next().vertex();

                if (visited[vid1] > 0 || visited[vid2] > 0 || visited[vid3] > 0)
                    continue;

                CGLA::Vec3d p1 = m.pos(vid1), p2 = m.pos(vid2),
                            p3 = m.pos(vid3);
                CGLA::Vec3d v1 = CGLA::cond_normalize(p2 - p1),
                            v2 = CGLA::cond_normalize(p3 - p2);

                // Compute the rotational angle and direction
                double angle = acos(dot(v1, v2));
                double sign =
                    sgn(dot(normal(m, w.opp().face()), cross(v1, v2)));

                // Check rotational angle and direction
                if (angle_threshold < angle || sign < 0.0) continue;

                // Add triangle, and pair it with the existing vertices
                // If new edge is too long, add new point and create triangle
                if (length(p3 - p1) > 8.0 * avg)
                    continue;
                else if (length(p3 - p1) > 4.0 * avg) {

                    CGLA::Vec3d p4 = 0.5 * (p1 + p3);
                    p4 = p2 + 0.5 * avg * CGLA::cond_normalize(p4 - p2);

                    HMesh::FaceID fid2 = m.add_face({p1, p2, p4});
                    HMesh::FaceID fid3 = m.add_face({p2, p3, p4});

                    HMesh::Walker w2 = m.walker(fid2), w3 = m.walker(fid3);

                    pair[w2.vertex()]        = pair[vid1];
                    pair[w2.next().vertex()] = pair[vid2];

                    pair[w3.vertex()]               = pair[vid2];
                    pair[w3.next().vertex()]        = pair[vid3];
                    pair[w3.next().next().vertex()] = pair_value++;

                } else {
                    HMesh::Walker w2 = m.walker(m.add_face({p1, p2, p3}));

                    pair[w2.vertex()]               = pair[vid1];
                    pair[w2.next().vertex()]        = pair[vid2];
                    pair[w2.next().next().vertex()] = pair[vid3];
                }
                visited[vid1] = 1, visited[vid2] = 1, visited[vid3] = 1;

                change = true;
            }

            // If triangles are added, then stitch them to the mesh
            if (!change) continue;

            stitch_mesh(m, pair);
            m.cleanup();
            HMesh::refine_edges(m, 2.0 * avg);
            HMesh::triangulate(m);
            m.cleanup();
        }
        m.cleanup();
    }
} // namespace Utilities

// End of File
// ========================================================================== //