#include <numeric>

#include <process.hpp>

#include <prim_gen.hpp>
#include <connect_by_topo.hpp>
#include <post_topo.hpp>

ISNP_API void build_implicit_network_by_blobtree(const s_settings&               settings,
                                                 const baked_blobtree_t&         tree,
                                                 stl_vector_mp<Eigen::Vector3d>& output_vertices,
                                                 stl_vector_mp<uint32_t>&        output_polygon_faces,
                                                 stl_vector_mp<uint32_t>&        output_vertex_counts_of_face)
{
    // load LUT for implicit arrangement
    load_lut();

    scene_bg_mesh_info_t                                               scene_bg_mesh_info{};
    stl_vector_mp<std::array<compact_bg_mesh_coord_t, 4>>              tetrahedrons{};
    flat_hash_map_mp<compact_bg_mesh_coord_t, compact_bg_mesh_coord_t> vertex_lexigraphical_adjacency{};
    stl_vector_mp<uint32_t>                                            tetrahedron_active_subface_start_index{};
    stl_vector_mp<uint32_t>                                            active_subface_indices{};
    stl_vector_mp<arrangement_t>                                       tetrahedron_arrangements{};
    flat_hash_map_mp<compact_bg_mesh_coord_t, stl_vector_mp<uint32_t>> zero_vertex_to_incident_tet_mapping{};
    stl_vector_mp<Eigen::Vector3d>                                     iso_pts{};
    stl_vector_mp<iso_vertex_t>                                        iso_verts{};
    stl_vector_mp<polygon_face_t>                                      iso_faces{};
    // primitive generation
    {
        flat_hash_map_mp<compact_bg_mesh_coord_t, bg_mesh_vert_info_t> vertex_infos{};
        {
            btree_map_mp<compact_bg_mesh_coord_t, stl_vector_mp<uint32_t>>                    vertex_to_tet_mapping{};
            flat_hash_map_mp<compact_bg_mesh_coord_t, stl_vector_mp<compact_bg_mesh_coord_t>> reverse_vertex_adjacency{};
            {
                stl_vector_mp<grid_region> vertex_exist_regions{};
                extract_vertex_infos(settings, tree, scene_bg_mesh_info, vertex_infos, vertex_exist_regions);
                build_tetrahedron_and_adjacency(scene_bg_mesh_info,
                                                vertex_exist_regions,
                                                tetrahedrons,
                                                vertex_lexigraphical_adjacency,
                                                reverse_vertex_adjacency,
                                                vertex_to_tet_mapping);
            }
            filter_tet_by_subface(vertex_to_tet_mapping,
                                  vertex_infos,
                                  tetrahedrons,
                                  vertex_lexigraphical_adjacency,
                                  reverse_vertex_adjacency,
                                  tetrahedron_active_subface_start_index,
                                  active_subface_indices,
                                  zero_vertex_to_incident_tet_mapping);
        }

        auto tet_active_subface_counts = compute_implicit_arrangements(vertex_infos,
                                                                       tetrahedrons,
                                                                       tetrahedron_active_subface_start_index,
                                                                       active_subface_indices,
                                                                       tetrahedron_arrangements);
        extract_iso_mesh(tet_active_subface_counts,
                         tetrahedron_arrangements,
                         scene_bg_mesh_info,
                         tetrahedrons,
                         tetrahedron_active_subface_start_index,
                         active_subface_indices,
                         vertex_infos,
                         iso_pts,
                         iso_verts,
                         iso_faces);
    }
    // connect components by topology
    // TODO: split headers of chains out, because it should not change during the chain
    stl_vector_mp<iso_edge_t> chain_representative_headers{};
    stl_vector_mp<uint32_t>   patch_of_face(iso_faces.size());
    stl_vector_mp<uint32_t>   shell_of_half_patch{};
    stl_vector_mp<uint32_t>   component_of_patch{};
    flat_index_group_t        patches{};
    flat_index_group_t        chains{};
    flat_index_group_t        shells{};
    flat_index_group_t        components{};
    flat_index_group_t        arrangement_cells{};
    stl_vector_mp<uint32_t>   shell_to_cell{};
    {
        {
            stl_vector_mp<stl_vector_mp<uint32_t>> edges_of_iso_face{};
            stl_vector_mp<iso_edge_t>              iso_edges{};
            compute_patch_edges(iso_faces, edges_of_iso_face, iso_edges);
            compute_patches(edges_of_iso_face, iso_edges, iso_faces, patches, patch_of_face);
            compute_chains(iso_verts.size(), iso_edges, chains, chain_representative_headers);
        }
        stl_vector_mp<stl_vector_mp<uint32_t>> half_patch_adj_list(2 * chains.size());
        chains.group_foreach([&](uint32_t chain_index, uint32_t _, uint32_t __) {
            compute_patch_order(tetrahedrons,
                                chain_representative_headers[chain_index],
                                iso_verts,
                                iso_faces,
                                tetrahedron_arrangements,
                                tetrahedron_active_subface_start_index,
                                active_subface_indices,
                                zero_vertex_to_incident_tet_mapping,
                                patch_of_face,
                                half_patch_adj_list);
        });
        compute_shells_and_components(half_patch_adj_list, shells, shell_of_half_patch, components, component_of_patch);

        if (components.size() == 1) // no nesting problem, each shell is an arrangement cell
        {
            arrangement_cells.index_group.resize(shells.size(), 0u);
            arrangement_cells.start_indices.resize(shells.size() + 1);
            std::iota(arrangement_cells.start_indices.begin(), arrangement_cells.start_indices.end(), 0u);
        } else {
            {
                stl_vector_mp<std::pair<uint32_t, uint32_t>> shell_links{};
                topo_ray_shooting(scene_bg_mesh_info,
                                  tetrahedrons,
                                  vertex_lexigraphical_adjacency,
                                  tetrahedron_arrangements,
                                  iso_verts,
                                  iso_faces,
                                  patches,
                                  patch_of_face,
                                  shells,
                                  shell_of_half_patch,
                                  components,
                                  component_of_patch,
                                  shell_links);
                compute_arrangement_cells(shells.size(), shell_links, arrangement_cells);
            }
        }

        shell_to_cell.resize(arrangement_cells.size());
        arrangement_cells.foreach ([&shell_to_cell](uint32_t cell_index, uint32_t _, uint32_t shell_index) {
            shell_to_cell[shell_index] = cell_index;
        });
    }

    // post process
    {
        dynamic_bitset_mp<> active_cell_label{};
        {
            stl_vector_mp<dynamic_bitset_mp<>> cell_primitive_signs(tree.primitives.size(),
                                                                    dynamic_bitset_mp<>(arrangement_cells.size()));
            {
                stl_vector_mp<dynamic_bitset_mp<>> cell_subface_signs(tree.subfaces.size(),
                                                                      dynamic_bitset_mp<>(arrangement_cells.size()));
                propagate_subface_labels(tree,
                                         iso_faces,
                                         patches,
                                         arrangement_cells,
                                         shell_of_half_patch,
                                         shells,
                                         shell_to_cell,
                                         cell_subface_signs);
                transform_subface_to_primitive_labels(tree, cell_subface_signs, cell_primitive_signs);
            }
            active_cell_label = filter_cells_by_boolean(tree, cell_primitive_signs);
        }
        filter_polygon_faces(iso_faces,
                             patches,
                             arrangement_cells,
                             shell_of_half_patch,
                             shells,
                             shell_to_cell,
                             active_cell_label,
                             output_polygon_faces,
                             output_vertex_counts_of_face);
        filter_active_vertices(output_vertices, output_polygon_faces);
    }
}