diff --git a/network_process/include/connect_by_topo/topology_ray_shooting.hpp b/network_process/include/connect_by_topo/topology_ray_shooting.hpp
index e60a59c..6f7e9ad 100644
--- a/network_process/include/connect_by_topo/topology_ray_shooting.hpp
+++ b/network_process/include/connect_by_topo/topology_ray_shooting.hpp
@@ -11,6 +11,7 @@ struct face_with_orient_t;
 // topological ray shooting for implicit arrangement
 void topo_ray_shooting(const scene_bg_mesh_info_t                   &scene_bg_mesh_info,
                        const stl_vector_mp<std::array<uint32_t, 4>> &tetrahedrons,
+                       const flat_hash_map_mp<uint32_t, uint32_t>   &vertex_lexigraphical_adjacency,
                        const stl_vector_mp<arrangement_t>           &tetrahedron_arrangements,
                        const stl_vector_mp<iso_vertex_t>            &iso_verts,
                        const stl_vector_mp<polygon_face_t>          &iso_faces,
@@ -22,12 +23,6 @@ void topo_ray_shooting(const scene_bg_mesh_info_t                   &scene_bg_me
                        const stl_vector_mp<uint32_t>                &component_of_patch,
                        stl_vector_mp<std::pair<uint32_t, uint32_t>> &shell_links);
 
-// Given tet mesh,
-// build the map: v-->v_next, where v_next has lower order than v
-void build_next_vert(const scene_bg_mesh_info_t                   &scene_bg_mesh_info,
-                     const stl_vector_mp<std::array<uint32_t, 4>> &tetrahedrons,
-                     flat_hash_map_mp<uint32_t, uint32_t>         &next_vert);
-
 // find extremal edge for each component
 void find_extremal_edges(const scene_bg_mesh_info_t                                       &scene_bg_mesh_info,
                          const stl_vector_mp<iso_vertex_t>                                &iso_verts,
diff --git a/network_process/include/prim_gen.hpp b/network_process/include/prim_gen.hpp
index 7606879..8901fd1 100644
--- a/network_process/include/prim_gen.hpp
+++ b/network_process/include/prim_gen.hpp
@@ -13,13 +13,17 @@ void                    extract_vertex_infos(const s_settings&           setting
                                              const baked_blobtree_t&     tree,
                                              scene_bg_mesh_info_t&       scene_bg_mesh_info,
                                              Eigen::Ref<Eigen::MatrixXd> vertex_infos);
-void                    pair_tetrahedron(const scene_bg_mesh_info_t&                          scene_bg_mesh_info,
-                                         stl_vector_mp<std::array<uint32_t, 4>>&              tetrahedrons,
-                                         flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& vertex_to_tet_mapping);
-void                    filter_tet_by_subface(const flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& vertex_to_tet_mapping,
+void                    build_tetrahedron_and_adjacency(const scene_bg_mesh_info_t&                          scene_bg_mesh_info,
+                                                        stl_vector_mp<std::array<uint32_t, 4>>&              tetrahedrons,
+                                                        flat_hash_map_mp<uint32_t, uint32_t>&                vertex_lexigraphical_adjacency,
+                                                        flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& reverse_vertex_adjacency,
+                                                        btree_map_mp<uint32_t, stl_vector_mp<uint32_t>>&     vertex_to_tet_mapping);
+void                    filter_tet_by_subface(const btree_map_mp<uint32_t, stl_vector_mp<uint32_t>>&     vertex_to_tet_mapping,
                                               Eigen::Ref<Eigen::MatrixXd>                                vertex_infos,
                                               flat_hash_map_mp<uint32_t, uint32_t>&                      vertex_indices_mapping,
                                               stl_vector_mp<std::array<uint32_t, 4>>&                    tetrahedrons,
+                                              flat_hash_map_mp<uint32_t, uint32_t>&                      vertex_lexigraphical_adjacency,
+                                              const flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& reverse_vertex_adjacency,
                                               stl_vector_mp<uint32_t>&                                   tetrahedron_active_subface_start_index,
                                               stl_vector_mp<uint32_t>&                                   active_subface_indices,
                                               flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>&       zero_vertex_to_incident_tet_mapping);
diff --git a/network_process/interface/fwd_types.hpp b/network_process/interface/fwd_types.hpp
index b899491..dbb6be1 100644
--- a/network_process/interface/fwd_types.hpp
+++ b/network_process/interface/fwd_types.hpp
@@ -2,6 +2,7 @@
 
 #include <container/dynamic_bitset.hpp>
 #include <container/hashmap.hpp>
+#include <container/btree.hpp>
 #include <container/stl_alias.hpp>
 #include <math/math_defs.hpp>
 
@@ -22,9 +23,10 @@ struct scene_bg_mesh_info_t {
     Eigen::Vector<uint64_t, 3> vert_resolution{};
 };
 
+// x as the most significant bit, y as the second most significant bit, z as the least significant bit
 static inline uint32_t hash_vert_pos(size_t x, size_t y, size_t z, const scene_bg_mesh_info_t& info)
 {
-    return static_cast<uint32_t>(x + info.vert_resolution.x() * (y + info.vert_resolution.y() * z));
+    return static_cast<uint32_t>(z + info.vert_resolution.z() * (y + info.vert_resolution.y() * x));
 }
 
 static inline Eigen::Vector3d get_grid_pos(size_t x, size_t y, size_t z, const scene_bg_mesh_info_t& info)
@@ -39,17 +41,17 @@ static inline Eigen::Vector3d get_grid_pos(const Eigen::Vector3d& pos, const sce
 
 static inline Eigen::Vector3d get_vert_pos(uint32_t hashed_value, const scene_bg_mesh_info_t& info)
 {
-    size_t x = static_cast<size_t>(hashed_value % (info.vert_resolution.x()));
-    size_t y = static_cast<size_t>((hashed_value / info.vert_resolution.x()) % info.vert_resolution.y());
-    size_t z = static_cast<size_t>(hashed_value / (info.vert_resolution.x() * info.vert_resolution.y()));
+    size_t z = static_cast<size_t>(hashed_value % (info.vert_resolution.z()));
+    size_t y = static_cast<size_t>((hashed_value / info.vert_resolution.z()) % info.vert_resolution.y());
+    size_t x = static_cast<size_t>(hashed_value / (info.vert_resolution.z() * info.vert_resolution.y()));
     return get_grid_pos(x, y, z, info);
 }
 
 static inline std::array<size_t, 3> get_grid_pos(uint32_t hashed_value, const scene_bg_mesh_info_t& info)
 {
-    size_t x = static_cast<size_t>(hashed_value % (info.vert_resolution.x()));
-    size_t y = static_cast<size_t>((hashed_value / info.vert_resolution.x()) % info.vert_resolution.y());
-    size_t z = static_cast<size_t>(hashed_value / (info.vert_resolution.x() * info.vert_resolution.y()));
+    size_t z = static_cast<size_t>(hashed_value % (info.vert_resolution.z()));
+    size_t y = static_cast<size_t>((hashed_value / info.vert_resolution.z()) % info.vert_resolution.y());
+    size_t x = static_cast<size_t>(hashed_value / (info.vert_resolution.z() * info.vert_resolution.y()));
     return {x, y, z};
 }
 
diff --git a/network_process/src/connect_by_topo/topology_ray_shooting.cpp b/network_process/src/connect_by_topo/topology_ray_shooting.cpp
index a6787a0..fcbada8 100644
--- a/network_process/src/connect_by_topo/topology_ray_shooting.cpp
+++ b/network_process/src/connect_by_topo/topology_ray_shooting.cpp
@@ -70,16 +70,11 @@ struct equal_to<face_with_orient_t> {
 };
 } // namespace std
 
-// point (x,y,z): dictionary order
-inline bool point_xyz_less(const Eigen::Vector3d &p, const Eigen::Vector3d &q)
-{
-    return std::lexicographical_compare(p.begin(), p.end(), q.begin(), q.end());
-}
-
 // ============================================================================================
 
 void topo_ray_shooting(const scene_bg_mesh_info_t                   &scene_bg_mesh_info,
                        const stl_vector_mp<std::array<uint32_t, 4>> &tetrahedrons,
+                       const flat_hash_map_mp<uint32_t, uint32_t>   &vertex_lexigraphical_adjacency,
                        const stl_vector_mp<arrangement_t>           &tetrahedron_arrangements,
                        const stl_vector_mp<iso_vertex_t>            &iso_verts,
                        const stl_vector_mp<polygon_face_t>          &iso_faces,
@@ -91,10 +86,6 @@ void topo_ray_shooting(const scene_bg_mesh_info_t                   &scene_bg_me
                        const stl_vector_mp<uint32_t>                &component_of_patch,
                        stl_vector_mp<std::pair<uint32_t, uint32_t>> &shell_links)
 {
-    // map: tet vert index --> index of next vert (with smaller (x,y,z))
-    flat_hash_map_mp<uint32_t, uint32_t> next_vert{};
-    build_next_vert(scene_bg_mesh_info, tetrahedrons, next_vert);
-
     // find extremal edge for each component
     // extremal edge of component i is stored at position [2*i], [2*i+1]
     stl_vector_mp<uint32_t>                                          extremal_edge_of_component{};
@@ -110,7 +101,7 @@ void topo_ray_shooting(const scene_bg_mesh_info_t                   &scene_bg_me
                         patches,
                         components,
                         component_of_patch,
-                        next_vert,
+                        vertex_lexigraphical_adjacency,
                         extremal_edge_of_component,
                         iso_vert_on_v_v_next,
                         iso_face_id_of_tet_face,
@@ -178,18 +169,16 @@ void topo_ray_shooting(const scene_bg_mesh_info_t                   &scene_bg_me
                                                                    : shell_of_half_patch[2 * patch_of_face[iso_fId] + 1];
                 // follow the ray till another iso-vertex or the sink
                 auto       v_curr      = extreme_v2;
-                // while (next_vert[v_curr] != invalid_index && iso_vert_on_v_v_next[v_curr] == invalid_index) {
-                while (next_vert.find(v_curr) != next_vert.end()
+                while (vertex_lexigraphical_adjacency.find(v_curr) != vertex_lexigraphical_adjacency.end()
                        && iso_vert_on_v_v_next.find(v_curr) == iso_vert_on_v_v_next.end()) {
-                    v_curr = next_vert[v_curr];
+                    v_curr = vertex_lexigraphical_adjacency.at(v_curr);
                 }
-                // if (iso_vert_on_v_v_next[v_curr] != invalid_index) {
                 if (iso_vert_on_v_v_next.find(v_curr) != iso_vert_on_v_v_next.end()) {
                     // reached iso-vert at end of the ray
                     const auto  iso_vId_end         = iso_vert_on_v_v_next[v_curr];
                     const auto  end_tetId           = iso_verts[iso_vId_end].header.volume_index;
                     const auto &end_tet_arrangement = tetrahedron_arrangements[end_tetId];
-                    auto        v_next              = next_vert[v_curr];
+                    auto        v_next              = vertex_lexigraphical_adjacency.at(v_curr);
                     // find local vertex indices in the end tetrahedron
                     for (uint32_t j = 0; j < 4; ++j) {
                         if (tetrahedrons[end_tetId][j] == v_curr) {
@@ -221,24 +210,6 @@ void topo_ray_shooting(const scene_bg_mesh_info_t                   &scene_bg_me
     }
 }
 
-void build_next_vert(const scene_bg_mesh_info_t                   &scene_bg_mesh_info,
-                     const stl_vector_mp<std::array<uint32_t, 4>> &tetrahedrons,
-                     flat_hash_map_mp<uint32_t, uint32_t>         &next_vert)
-{
-    next_vert.reserve(tetrahedrons.size() * 2);
-
-    for (const auto &tet : tetrahedrons) {
-        // find the smallest vertex of tet
-        auto min_iter = std::min_element(tet.begin(), tet.end(), [&](const uint32_t &lhs, const uint32_t &rhs) {
-            return point_xyz_less(get_vert_pos(lhs, scene_bg_mesh_info), get_vert_pos(rhs, scene_bg_mesh_info));
-        });
-        //
-        for (auto iter = tet.begin(); iter != tet.end(); ++iter) {
-            if (iter != min_iter) next_vert[*iter] = *min_iter;
-        }
-    }
-}
-
 void find_extremal_edges(const scene_bg_mesh_info_t                                       &scene_bg_mesh_info,
                          const stl_vector_mp<iso_vertex_t>                                &iso_verts,
                          const stl_vector_mp<polygon_face_t>                              &iso_faces,
@@ -277,12 +248,8 @@ void find_extremal_edges(const scene_bg_mesh_info_t
                                 u2 = v2;
                             } else {
                                 if (v2 == u2) {
-                                    if (point_xyz_less(get_vert_pos(v1, scene_bg_mesh_info),
-                                                       get_vert_pos(u1, scene_bg_mesh_info))) {
-                                        u1 = v1;
-                                    }
-                                } else if (point_xyz_less(get_vert_pos(v2, scene_bg_mesh_info),
-                                                          get_vert_pos(u2, scene_bg_mesh_info))) {
+                                    if (v1 < u1) { u1 = v1; }
+                                } else if (v2 < u2) {
                                     u1 = v1;
                                     u2 = v2;
                                 }
@@ -300,12 +267,8 @@ void find_extremal_edges(const scene_bg_mesh_info_t
                                 u2 = v1;
                             } else {
                                 if (v1 == u2) {
-                                    if (point_xyz_less(get_vert_pos(v2, scene_bg_mesh_info),
-                                                       get_vert_pos(u1, scene_bg_mesh_info))) {
-                                        u1 = v2;
-                                    }
-                                } else if (point_xyz_less(get_vert_pos(v1, scene_bg_mesh_info),
-                                                          get_vert_pos(u2, scene_bg_mesh_info))) {
+                                    if (v2 < u1) { u1 = v2; }
+                                } else if (v1 < u2) {
                                     u1 = v2;
                                     u2 = v1;
                                 }
diff --git a/network_process/src/prim_gen/pair_tetrahedron.cpp b/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp
similarity index 72%
rename from network_process/src/prim_gen/pair_tetrahedron.cpp
rename to network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp
index f17ed9a..25247b4 100644
--- a/network_process/src/prim_gen/pair_tetrahedron.cpp
+++ b/network_process/src/prim_gen/build_tetrahedron_and_adjacency.cpp
@@ -2,9 +2,24 @@
 
 static constexpr uint64_t intertwined_01_bit_block = 0x5555555555555555;
 
-void pair_tetrahedron(const scene_bg_mesh_info_t&                          scene_bg_mesh_info,
-                      stl_vector_mp<std::array<uint32_t, 4>>&              tetrahedrons,
-                      flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& vertex_to_tet_mapping)
+static void insert_or_compare_vertex_adjacency(flat_hash_map_mp<uint32_t, uint32_t>& vertex_lexigraphical_adjacency,
+                                               flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& reverse_vertex_adjacency,
+                                               std::array<uint32_t, 3>                              src,
+                                               uint32_t                                             dst)
+{
+    reverse_vertex_adjacency[dst].reserve(reverse_vertex_adjacency[dst].size() + 3);
+    for (const auto& s : src) {
+        auto [iter, is_new] = vertex_lexigraphical_adjacency.try_emplace(s, dst);
+        if (!is_new) iter->second = std::min(iter->second, dst);
+        reverse_vertex_adjacency[dst].emplace_back(s);
+    }
+}
+
+void build_tetrahedron_and_adjacency(const scene_bg_mesh_info_t&                          scene_bg_mesh_info,
+                                     stl_vector_mp<std::array<uint32_t, 4>>&              tetrahedrons,
+                                     flat_hash_map_mp<uint32_t, uint32_t>&                vertex_lexigraphical_adjacency,
+                                     flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& reverse_vertex_adjacency,
+                                     btree_map_mp<uint32_t, stl_vector_mp<uint32_t>>&     vertex_to_tet_mapping)
 {
     auto max_pos = scene_bg_mesh_info.vert_resolution.template cast<size_t>() - Eigen::Vector<size_t, 3>::Ones();
     tetrahedrons.reserve((max_pos.array() - 1).prod() * 5);
@@ -27,12 +42,52 @@ void pair_tetrahedron(const scene_bg_mesh_info_t&                          scene
                     tetrahedrons.emplace_back(std::array{v1, v3, v0, v4});
                     tetrahedrons.emplace_back(std::array{v3, v1, v2, v6});
                     tetrahedrons.emplace_back(std::array{v4, v1, v6, v5});
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v6, v1, v3},
+                                                       v4);
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v6, v3, v7},
+                                                       v4);
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v1, v3, v4},
+                                                       v0);
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v1, v2, v6},
+                                                       v3);
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v1, v6, v5},
+                                                       v4);
                 } else {
                     tetrahedrons.emplace_back(std::array{v7, v0, v2, v5});
                     tetrahedrons.emplace_back(std::array{v2, v3, v0, v7});
                     tetrahedrons.emplace_back(std::array{v5, v7, v0, v4});
                     tetrahedrons.emplace_back(std::array{v7, v2, v6, v5});
                     tetrahedrons.emplace_back(std::array{v0, v1, v2, v5});
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v7, v2, v5},
+                                                       v0);
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v2, v3, v7},
+                                                       v0);
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v5, v7, v4},
+                                                       v0);
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v2, v6, v5},
+                                                       v7);
+                    insert_or_compare_vertex_adjacency(vertex_lexigraphical_adjacency,
+                                                       reverse_vertex_adjacency,
+                                                       {v1, v2, v5},
+                                                       v0);
                 }
             }
         }
diff --git a/network_process/src/prim_gen/extract_mesh.cpp b/network_process/src/prim_gen/extract_mesh.cpp
index 004a523..39bd83e 100644
--- a/network_process/src/prim_gen/extract_mesh.cpp
+++ b/network_process/src/prim_gen/extract_mesh.cpp
@@ -7,7 +7,7 @@ void extract_iso_mesh(const std::array<uint32_t, 3>&                tet_active_s
                       const stl_vector_mp<std::array<uint32_t, 4>>& tetrahedrons,
                       const stl_vector_mp<uint32_t>&                tetrahedron_active_subface_start_index,
                       const stl_vector_mp<uint32_t>&                active_subface_indices,
-                      const Eigen::MatrixXd&                        vertex_infos,
+                      const Eigen::Ref<const Eigen::MatrixXd>       vertex_infos,
                       const flat_hash_map_mp<uint32_t, uint32_t>&   vertex_indices_mapping,
                       stl_vector_mp<Eigen::Vector3d>&               iso_pts,
                       stl_vector_mp<iso_vertex_t>&                  iso_verts,
diff --git a/network_process/src/prim_gen/filter_tet_by_subface.cpp b/network_process/src/prim_gen/filter_tet_by_subface.cpp
index 2a83adf..7f44c3f 100644
--- a/network_process/src/prim_gen/filter_tet_by_subface.cpp
+++ b/network_process/src/prim_gen/filter_tet_by_subface.cpp
@@ -15,10 +15,12 @@
 //     }
 // };
 
-void filter_tet_by_subface(const flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& vertex_to_tet_mapping,
+void filter_tet_by_subface(const btree_map_mp<uint32_t, stl_vector_mp<uint32_t>>&     vertex_to_tet_mapping,
                            Eigen::Ref<Eigen::MatrixXd>                                vertex_infos,
                            flat_hash_map_mp<uint32_t, uint32_t>&                      vertex_indices_mapping,
                            stl_vector_mp<std::array<uint32_t, 4>>&                    tetrahedrons,
+                           flat_hash_map_mp<uint32_t, uint32_t>&                      vertex_lexigraphical_adjacency,
+                           const flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>& reverse_vertex_adjacency,
                            stl_vector_mp<uint32_t>&                                   tetrahedron_active_subface_start_index,
                            stl_vector_mp<uint32_t>&                                   active_subface_indices,
                            flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>>&       zero_vertex_to_incident_tet_mapping)
@@ -88,16 +90,57 @@ void filter_tet_by_subface(const flat_hash_map_mp<uint32_t, stl_vector_mp<uint32
     filtered_active_vertices.shrink_to_fit();
     std::swap(filtered_tets, tetrahedrons);
 
+    // change vertex adjacency with filtered vertices, so that each vertex pair in adjacency map is connected directly
+    // without any intermediate vertices
     std::sort(filtered_active_vertices.begin(), filtered_active_vertices.end());
-    auto            end_iter     = std::unique(filtered_active_vertices.begin(), filtered_active_vertices.end());
-    auto            vertex_count = std::distance(filtered_active_vertices.begin(), end_iter);
-    Eigen::MatrixXd filtered_vert_infos(vertex_count, vertex_infos.cols());
-    auto            filtered_vert_infos_iter = filtered_vert_infos.rowwise().begin();
-    for (auto iter = filtered_active_vertices.begin(); iter != end_iter; ++iter) {
+    auto                     end_iter = std::unique(filtered_active_vertices.begin(), filtered_active_vertices.end());
+    auto                     iter     = filtered_active_vertices.begin();
+    stl_vector_mp<uint32_t*> prev_to_next_adjacency{};
+    uint32_t                 next_vertex_index{};
+    prev_to_next_adjacency.reserve(3);
+    for (auto i = 0; i < vertex_infos.rows(); ++i) {
+        if (*iter == i) { // the vertex is still active, skip it
+            ++iter;
+        } else {          // the vertex is filtered out, reconnect its prev to next
+            prev_to_next_adjacency.clear();
+            // the vertex at upper bound do not have reversed adjacency
+            if (auto iter = reverse_vertex_adjacency.find(i); iter != reverse_vertex_adjacency.end()) {
+                for (const auto& prev : iter->second) {
+                    prev_to_next_adjacency.emplace_back(&vertex_lexigraphical_adjacency[prev]);
+                }
+            }
+            // also, the vertex at lower bound may not have adjacency
+            if (auto iter = vertex_lexigraphical_adjacency.find(i); iter != vertex_lexigraphical_adjacency.end()) {
+                next_vertex_index = iter->second;
+            } else {
+                next_vertex_index = std::numeric_limits<uint32_t>::max();
+            }
+            // reconnect all previous vertices to next vertex
+            for (auto prev : prev_to_next_adjacency) { *prev = next_vertex_index; }
+        }
+    }
+    iter = filtered_active_vertices.begin();
+    // remove filtered vertices from vertex_lexigraphical_adjacency
+    for (auto i = 0; i < vertex_infos.rows(); ++i) {
+        if (*iter == i) {
+            ++iter;
+        } else {
+            vertex_lexigraphical_adjacency.erase(i);
+        }
+    }
+
+    auto                                 vertex_count = std::distance(filtered_active_vertices.begin(), end_iter);
+    Eigen::MatrixXd                      filtered_vert_infos(vertex_count, vertex_infos.cols());
+    flat_hash_map_mp<uint32_t, uint32_t> filtered_vertex_lexigraphical_adjacency{};
+    auto                                 filtered_vert_infos_iter = filtered_vert_infos.rowwise().begin();
+    iter                                                          = filtered_active_vertices.begin();
+    for (; iter != end_iter; ++iter) {
         const auto vertex_index              = *iter;
         *filtered_vert_infos_iter            = vertex_infos.row(vertex_index);
         vertex_indices_mapping[vertex_index] = std::distance(filtered_active_vertices.begin(), iter);
+        filtered_vertex_lexigraphical_adjacency.emplace(vertex_index, vertex_indices_mapping.at(vertex_index));
         filtered_vert_infos_iter++;
     }
-    vertex_infos = filtered_vert_infos;
+    vertex_infos                   = std::move(filtered_vert_infos);
+    vertex_lexigraphical_adjacency = std::move(filtered_vertex_lexigraphical_adjacency);
 }
\ No newline at end of file
diff --git a/network_process/src/process.cpp b/network_process/src/process.cpp
index e6a8330..8a8bc33 100644
--- a/network_process/src/process.cpp
+++ b/network_process/src/process.cpp
@@ -15,6 +15,7 @@ ISNP_API void build_implicit_network_by_blobtree(const s_settings&
 
     scene_bg_mesh_info_t                                scene_bg_mesh_info{};
     stl_vector_mp<std::array<uint32_t, 4>>              tetrahedrons{};
+    flat_hash_map_mp<uint32_t, uint32_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{};
@@ -27,16 +28,23 @@ ISNP_API void build_implicit_network_by_blobtree(const s_settings&
         Eigen::MatrixXd                      vertex_infos{};
         flat_hash_map_mp<uint32_t, uint32_t> vertex_indices_mapping{};
         {
-            flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>> vertex_to_tet_mapping{};
+            btree_map_mp<uint32_t, stl_vector_mp<uint32_t>>     vertex_to_tet_mapping{};
+            flat_hash_map_mp<uint32_t, stl_vector_mp<uint32_t>> reverse_vertex_adjacency{};
             {
                 stl_vector_mp<grid_region> vertex_exist_regions{};
                 extract_vertex_infos(settings, tree, scene_bg_mesh_info, vertex_infos);
-                pair_tetrahedron(scene_bg_mesh_info, tetrahedrons, vertex_to_tet_mapping);
+                build_tetrahedron_and_adjacency(scene_bg_mesh_info,
+                                                tetrahedrons,
+                                                vertex_lexigraphical_adjacency,
+                                                reverse_vertex_adjacency,
+                                                vertex_to_tet_mapping);
             }
             filter_tet_by_subface(vertex_to_tet_mapping,
                                   vertex_infos,
                                   vertex_indices_mapping,
                                   tetrahedrons,
+                                  vertex_lexigraphical_adjacency,
+                                  reverse_vertex_adjacency,
                                   tetrahedron_active_subface_start_index,
                                   active_subface_indices,
                                   zero_vertex_to_incident_tet_mapping);
@@ -68,6 +76,8 @@ ISNP_API void build_implicit_network_by_blobtree(const s_settings&
     stl_vector_mp<stl_list_mp<iso_edge_t>> chains{};
     stl_vector_mp<stl_vector_mp<uint32_t>> shells{};
     stl_vector_mp<stl_vector_mp<uint32_t>> components{};
+    stl_vector_mp<stl_vector_mp<uint32_t>> arrangement_cells{};
+    stl_vector_mp<uint32_t>                shell_to_cell{};
     {
         {
             stl_vector_mp<stl_vector_mp<uint32_t>> edges_of_iso_face{};
@@ -89,43 +99,36 @@ ISNP_API void build_implicit_network_by_blobtree(const s_settings&
                                 patch_of_face,
                                 half_patch_adj_list);
         compute_shells_and_components(half_patch_adj_list, shells, shell_of_half_patch, components, component_of_patch);
-    }
-    // mixed topology connection and post process
-    if (components.size() == 1) // no nesting problem, each shell is an arrangement cell
-    {
-        output_vertices = std::move(iso_pts);
-        for (const auto& half_patch : shells[1]) {
-            for (const auto& face_id : patches[half_patch / 2]) {
-                const auto& face = iso_faces[face_id];
-                output_polygon_faces.insert(output_polygon_faces.end(),
-                                            std::make_move_iterator(face.vertex_indices.begin()),
-                                            std::make_move_iterator(face.vertex_indices.end()));
-                output_vertex_counts_of_face.emplace_back(face.vertex_indices.size());
-            }
-        }
-    } else {
-        stl_vector_mp<stl_vector_mp<uint32_t>> arrangement_cells{};
-        {
-            stl_vector_mp<std::pair<uint32_t, uint32_t>> shell_links{};
-            topo_ray_shooting(scene_bg_mesh_info,
-                              tetrahedrons,
-                              tetrahedron_arrangements,
-                              iso_verts,
-                              iso_faces,
-                              patches,
-                              patch_of_face,
-                              shells,
-                              shell_of_half_patch,
-                              components,
-                              component_of_patch,
-                              shell_links);
-            compute_arrangement_cells(static_cast<uint32_t>(shells.size()), shell_links, arrangement_cells);
-        }
+
+        if (components.size() == 1) // no nesting problem, each shell is an arrangement cell
         {
-            stl_vector_mp<uint32_t> shell_to_cell(shells.size());
+            arrangement_cells.reserve(shells.size());
+            for (uint32_t i = 0; i < shells.size(); ++i) { arrangement_cells.emplace_back(stl_vector_mp<uint32_t>{i}); }
+        } 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(static_cast<uint32_t>(shells.size()), shell_links, arrangement_cells);
+            }
+            shell_to_cell.resize(shells.size());
             for (uint32_t i = 0; i < arrangement_cells.size(); i++) {
                 for (auto shell : arrangement_cells[i]) shell_to_cell[shell] = i;
             }
+        }
+        // post process
+        {
             dynamic_bitset_mp<> active_cell_label{};
             {
                 stl_vector_mp<dynamic_bitset_mp<>> cell_primitive_signs(tree.primitives.size(),