// ========================================================================== //
// Class inclusion
#include "StreamSurfaces/StreamSurfaces.hpp"

// Local libraries
#include "Frames/FrameField.hpp"
#include "Utilities/MeshDistance.h"
#include "Pointcloud/Pointcloud.hpp"
#include "Utilities/math_operations.hpp"
#include "Utilities/mesh_operations.hpp"
#include "Utilities/string_operations.hpp"

// External libraries
#include <GEL/CGLA/CGLA.h>
#include <GEL/Geometry/RGrid.h>
#include <GEL/Geometry/XForm.h>
#include <GEL/HMesh/HMesh.h>
#include <GEL/HMesh/obj_load.h>
#include <GEL/HMesh/polygonize.h>

// Standard libraries
#include <filesystem>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>

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

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

// ========================================================================== //
// Wrap

HMesh::Manifold StreamSurfaces::isocontour_hull() {

    Utilities::debug_info(__FILE__, __LINE__, "Iso contouring", debug);

    // Grid resolution
    double      scale = 1.0;
    CGLA::Vec3i dims  = CGLA::Vec3i(scale * bbox.diagonal()) + CGLA::Vec3i(1);

    // Setup the xform and grid classes
    Geometry::XForm        xform(bbox.min_corner, bbox.max_corner, dims, 0.0);
    Geometry::RGrid<float> grid(dims, -1.0);

    // Populate the grid
    for (CGLA::Vec3i pi : Geometry::Range3D(dims)) {
        CGLA::Vec3d p = xform.inverse(CGLA::Vec3d(pi));

        p = CGLA::v_min(CGLA::v_max(p, bbox.min_corner), bbox.max_corner);
        std::vector<double> layer   = frame->interpolate_layer(p);
        double              density = frame->interpolate_density(p);

        double value_l = Utilities::max(layer) - threshold_void;
        double value_d = density - threshold_solid;

        double value = std::max(value_l, value_d);
        grid.store(pi, value);
    }

    // Compute the Iso-Contour
    HMesh::Manifold mesh;
    HMesh::volume_polygonize(xform, grid, mesh, 0.0);

    Utilities::debug_info(__FILE__, __LINE__, "Iso contouring Complete", debug);

    return mesh;
}

void StreamSurfaces::prepare_fixed() {

    if (!include_hull) return;
    Utilities::debug_info(__FILE__, __LINE__, "Preparing fixed", debug);

    HMesh::Manifold mesh = isocontour_hull();

    StreamSurface surf;
    for (HMesh::VertexID vid : mesh.vertices()) {
        CGLA::Vec3d         p     = mesh.pos(vid);
        CGLA::Vec3d         n     = HMesh::normal(mesh, vid);
        std::vector<double> layer = frame->interpolate_layer(p);
        double              l = std::max(Utilities::max(layer), threshold_void);

        // Add the mesh to the surfaces
        surf.pcl.push_back(p, n, l);
    }

    surf.mesh = mesh;
    fixed.push_back((int)surfaces.size());
    surfaces.push_back(surf);
}

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