ImplicitSurfaceNetwork/network_process/src/post_topo/patch_propagation.cpp

#include <stack>
#include <queue>

#include <post_topo.hpp>

void propagate_subface_labels(const baked_blobtree_t&                       tree,
                              const stl_vector_mp<polygon_face_t>&          faces,
                              const stl_vector_mp<stl_vector_mp<uint32_t>>& patches,
                              const stl_vector_mp<stl_vector_mp<uint32_t>>& arrangement_cells,
                              const stl_vector_mp<uint32_t>&                shell_of_half_patch,
                              const stl_vector_mp<stl_vector_mp<uint32_t>>& shells,
                              const stl_vector_mp<uint32_t>&                shell_to_cell,
                              stl_vector_mp<dynamic_bitset_mp<>>&           cell_subface_signs)
{
    const auto num_subface = tree.subfaces.size();

    stl_vector_mp<bool>                                   visited_cells(arrangement_cells.size(), false);
    stl_vector_mp<bool>                                   visited_subfaces(num_subface, false);
    stl_vector_mp<stl_vector_mp<uint32_t>>                cell_indices_of_inactive_subfaces(num_subface);
    std::queue<uint32_t>                                  Q{};
    flat_hash_map_mp<uint32_t, std::pair<uint32_t, bool>> cell_neighbors_map{};

    Q.emplace(0);
    while (!Q.empty()) {
        const auto cell_index = Q.front();
        const auto cell       = arrangement_cells[cell_index];
        Q.pop();

        if (!visited_cells[cell_index]) {
            visited_cells[cell_index] = true;
            cell_neighbors_map.clear();

            for (auto shell : cell) {
                for (auto half_patch : shells[shell]) {
                    auto subface_label = faces[patches[half_patch / 2][0]].subface_index;
                    // CAUTION: we assume that the sign is 1 when the surface is inside the sdf before
                    // but in blobtree we assume that the sign is 0 when the surface is inside the sdf
                    // so we need to flip the sign here
                    // bool sign           = (half_patch % 2 == 0) ? 1 : 0;
                    // auto oppose_cell    = shell_to_cell[shell_of_half_patch[sign ? (half_patch + 1) : (half_patch - 1)]];
                    bool sign          = (half_patch % 2 == 0) ? 0 : 1;
                    auto oppose_cell   = shell_to_cell[shell_of_half_patch[!sign ? (half_patch + 1) : (half_patch - 1)]];
                    cell_neighbors_map[oppose_cell] = std::make_pair(subface_label, sign);

#ifndef RELEASE_BRANCH
                    if (visited_subfaces[subface_label] != false && cell_subface_signs[subface_label][cell_index] != sign) {
                        throw std::runtime_error("ERROR: Inconsistent Cell Function Labels.");
                    }
#endif

                    cell_subface_signs[subface_label][cell_index] = sign;
                    visited_subfaces[subface_label]               = true;
                    // propagate the function signs to all previously inactive cells
                    if (cell_indices_of_inactive_subfaces[subface_label].size() > 0) {
                        for (const auto cell_index : cell_indices_of_inactive_subfaces[subface_label]) {
                            cell_subface_signs[subface_label][cell_index] = sign;
                        }
                        cell_indices_of_inactive_subfaces[subface_label].clear();
                    }
                }
            }

            // fetch inactive subface index
            for (uint32_t subface_index = 0; subface_index < num_subface; subface_index++) {
                if (visited_subfaces[subface_index] == false) {
                    cell_indices_of_inactive_subfaces[subface_index].emplace_back(cell_index);
                    for (const auto& [other_cell_index, _] : cell_neighbors_map) {
                        cell_indices_of_inactive_subfaces[subface_index].emplace_back(other_cell_index);
                    }
                }
            }

            // propagate to neighboring cells
            for (const auto& [other_cell_index, other_cell_func_label] : cell_neighbors_map) {
                if (!visited_cells[other_cell_index]) Q.emplace(other_cell_index);

                const auto& [func_index, sign] = other_cell_func_label;
                for (uint32_t subface_index = 0; subface_index < func_index; ++subface_index)
                    cell_subface_signs[subface_index][other_cell_index] = cell_subface_signs[subface_index][cell_index];
                // opposite cell has opposite sign on same patch
                cell_subface_signs[func_index][other_cell_index] = !sign;
                for (uint32_t subface_index = func_index + 1; subface_index < num_subface; ++subface_index)
                    cell_subface_signs[subface_index][other_cell_index] = cell_subface_signs[subface_index][cell_index];
            }
        }
    }

    // ASSUME: all subfaces should generate polygon mesh accordingly, so there should be no unknown subface signs
    //
    // // for all unvisited functions, they must totally contain or be contained by some surfaces
    // // so we can test their signs by testing whether the representative vertex is inside or outside the aabbs
    // for (uint32_t i = 0; i < num_subface; i++) {
    //     if (visited_subfaces[i] == false) {
    //         for (uint32_t j = 0; j < arrangement_cells.size(); ++j) {
    //             const auto& cell                  = arrangement_cells[j];
    //             const auto& representative_shell  = shells[cell[0]];
    //             const auto& representative_patch  = patches[representative_shell[0] / 2];
    //             const auto& representative_face   = faces[representative_patch[0]];
    //             const auto& representative_vertex = vertices[representative_face.vertex_indices[0]];

    //             cell_subface_signs[i][j] = get_primitive_node(i).aabb.contains(representative_vertex);
    //         }
    //     }
    // }
#ifndef RELEASE_BRANCH
    for (size_t i = 0; i < num_subface; ++i) {
        if (visited_subfaces[i] == false) { throw std::logic_error("ERROR: Still have sign-unknown subfaces."); }
    }
#endif
}

void transform_subface_to_primitive_labels(const baked_blobtree_t&                   tree,
                                           const stl_vector_mp<dynamic_bitset_mp<>>& cell_subface_signs,
                                           stl_vector_mp<dynamic_bitset_mp<>>&       cell_primitive_signs)
{   
    std::cout << "transform_subface_to_primitive_labels debug info: " << std::endl;
    for (size_t i = 0; i < cell_subface_signs.size(); ++i) {
        std::cout << "Subface " << i << " subface signs: " << std::endl;
        for (size_t j = 0; j < cell_subface_signs[i].size(); ++j) {
            std::cout << "  subface sign for cell " << j << ": " << cell_subface_signs[i][j] << std::endl;
        }
    }
    for (size_t j = 0; j < cell_subface_signs[0].size(); ++j) {
        std::cout << "Cell " << j << " subface signs: " ;
        for (size_t i = 0; i < cell_subface_signs.size(); ++i) {
            std::cout << " " << cell_subface_signs[i][j];
        }
        std::cout << std::endl;
    }
    stl_vector_mp<dynamic_bitset_mp<>> temp_subface_signs{};
    for (size_t i = 0; i < tree.primitives.size(); ++i) {
        auto&       cell_primitive_sign = cell_primitive_signs[i];
        const auto& primitive           = *tree.primitives[i].object_ptr;
        const auto& subface_indices     = tree.primitives[i].index_mapping;

        temp_subface_signs.clear();
        temp_subface_signs.reserve(subface_indices.size());
        for (const auto& subface_index : subface_indices) 
            temp_subface_signs.emplace_back(cell_subface_signs[subface_index]);
        // sign 0 in blobtree means inside, so we need to flip the sign
        cell_primitive_sign = primitive.judge_sign_by_subface_sign(temp_subface_signs);
        std::cout << "primitive " << i << " sign: " << cell_primitive_sign << std::endl;
        for (size_t j = 0; j < cell_primitive_sign.size(); ++j) {
            std::cout << "  subface sign for cell " << j << ": " << cell_primitive_sign[j] << std::endl;
        }
    }
}

dynamic_bitset_mp<> filter_cells_by_boolean(const baked_blobtree_t&             tree,
                                            stl_vector_mp<dynamic_bitset_mp<>>& cell_primitive_signs)
{
    struct compact_node_info {
        dynamic_bitset_mp<> cell_signs{};
        uint32_t            parent_index{};
    };

    std::stack<compact_node_info> stacked_nodes{};
    auto                          iter = tree.nodes.begin();
    compact_node_info             front_info{std::move(cell_primitive_signs[iter->primitive_index]), iter->parent_index};
    stacked_nodes.emplace(std::move(front_info));

    iter++;
    while (iter != tree.nodes.end()) {
        // each out iteration must start with leaf node, only 1 leaf node is absorbed in one iteration
        assert(iter->is_primitive_node());
        compact_node_info temp_info{std::move(cell_primitive_signs[iter->primitive_index]), iter->parent_index};
        iter++; // to parent or neighboring node
        while (!stacked_nodes.empty() && temp_info.parent_index == stacked_nodes.top().parent_index) {
            // do bool operation, util meet next primitive node.
            assert(iter->is_operation_node());

            const auto& other_cell_sign = stacked_nodes.top().cell_signs;
            switch (iter->get_operation()) {
                case internal::eNodeOperation::unionOp:        temp_info.cell_signs &= other_cell_sign; break;
                case internal::eNodeOperation::intersectionOp: temp_info.cell_signs |= other_cell_sign; break;
                case internal::eNodeOperation::differenceOp:
                    // stacked nodes are always left childs
                    temp_info.cell_signs.flip() |= other_cell_sign;
                    break;
                default: throw std::runtime_error("ERROR: baked blobtree with unknown type operation node"); break;
            }
            temp_info.parent_index = iter->parent_index;

            stacked_nodes.pop();
            iter++; // to parent or neighboring node
        }
        stacked_nodes.emplace(std::move(temp_info));
    }

    assert(stacked_nodes.size() == 1);
    assert(stacked_nodes.top().parent_index == 0xFFFFFFFF);

    // sign 0 in blobtree means inside, so we need to flip the sign
    return stacked_nodes.top().cell_signs.flip();
}