diff --git a/src/bonsai/bonsai/bim/module/model/polyline.py b/src/bonsai/bonsai/bim/module/model/polyline.py
index e33fe2cf5d..be491e756b 100644
--- a/src/bonsai/bonsai/bim/module/model/polyline.py
+++ b/src/bonsai/bonsai/bim/module/model/polyline.py
@@ -421,7 +421,7 @@ class PolylineOperator:
                     tool.Polyline.calculate_x_y_and_z(context, self.input_ui, self.tool_state)
 
                 tool.Blender.update_viewport()
-                return {"RUNNING_MODAL"}
+            return {"RUNNING_MODAL"}
 
     def set_offset(self, context: bpy.types.Context, relating_type: ifcopenshell.entity_instance) -> None:
         props = tool.Model.get_model_props()
@@ -461,6 +461,7 @@ class PolylineOperator:
         self.tool_state.axis_method = None
         self.tool_state.plane_method = None
         self.tool_state.mode = "Mouse"
+        tool.Raycast.clear_snap_objs()
         self.visible_objs = tool.Raycast.get_visible_objects(context)
         for obj in self.visible_objs:
             if bbox_2d := tool.Raycast.get_on_screen_2d_bounding_boxes(context, obj):
diff --git a/src/bonsai/bonsai/tool/raycast.py b/src/bonsai/bonsai/tool/raycast.py
index d98aaf26ab..0248dfeb56 100644
--- a/src/bonsai/bonsai/tool/raycast.py
+++ b/src/bonsai/bonsai/tool/raycast.py
@@ -16,6 +16,9 @@
 # You should have received a copy of the GNU General Public License
 # along with Bonsai.  If not, see <http://www.gnu.org/licenses/>.
 
+from __future__ import annotations
+
+import math
 from typing import Union
 
 import bmesh
@@ -41,6 +44,7 @@ class Raycast(bonsai.core.tool.Raycast):
         (0, -offset),
         (offset, -offset),
     )
+    snap_objs = []
 
     @classmethod
     def get_visible_objects(cls, context: bpy.types.Context):
@@ -232,6 +236,146 @@ class Raycast(bonsai.core.tool.Raycast):
         else:
             return None, None, None
 
+    @classmethod
+    def ray_cast_by_proximity_2d(
+        cls,
+        context: bpy.types.Context,
+        event: bpy.types.Event,
+        snap_obj,
+    ):
+        # TODO Clean unnecessary code
+
+        def divide_vector(start, end, n):
+            points = []
+            delta = (end - start) / n
+            for i in range(1, n):
+                point = start + i * delta
+                points.append(point)
+            return points
+
+        region = context.region
+        rv3d = context.region_data
+        mouse_pos = event.mouse_region_x, event.mouse_region_y
+        ray_origin, ray_target, ray_direction = cls.get_viewport_ray_data(context, event)
+        points = []
+
+        verts_2d = [
+            view3d_utils.location_3d_to_region_2d(region, rv3d, v) for v in snap_obj.verts_3d
+        ]  # Numpy version is worst in performance
+        verts_2d = [
+            view3d_utils.location_3d_to_region_2d(region, rv3d, v) for v in snap_obj.verts_3d
+        ]  # Numpy version is worst in performance
+        snap_threshold = 10.0
+
+        try:
+            loc = tool.Cad.region_2d_to_location_3d_np(region, rv3d, mouse_pos, ray_direction)
+        except:
+            loc = Vector((0, 0, 0))
+
+        for i, point in enumerate(verts_2d):
+            if not point:
+                break
+            distance = (Vector(mouse_pos) - point).length
+            if distance <= snap_threshold:
+                snap_point = {
+                    "object": snap_obj.obj,
+                    "type": "Vertex",
+                    "point": snap_obj.verts_3d[i],
+                    "distance": distance / 10,
+                }
+                points.append(snap_point)
+        intersected = snap_obj.raycast_boxes(
+            context, event, snap_obj.root, intersected=[], rays=(ray_origin, ray_direction)
+        )
+        edges = []
+        for it in intersected:
+            edges.extend(it.edges)
+        edges = set(edges)
+        count = 0
+        distances = {}
+        for e in edges:
+            idx = snap_obj.obj.data.edges[e].vertices
+
+            p0x, p0y = verts_2d[idx[0]][0], verts_2d[idx[0]][1]
+            p1x, p1y = verts_2d[idx[1]][0], verts_2d[idx[1]][1]
+            px, py = mouse_pos
+
+            # segment vector = p1 - p0
+            sx = p1x - p0x
+            sy = p1y - p0y
+
+            # seg length squared
+            seg_len_sq = sx * sx + sy * sy
+
+            if seg_len_sq == 0.0:
+                # degenerate segment: return distance to p0
+                dx = px - p0x
+                dy = py - p0y
+                dist = math.hypot(dx, dy)
+                return dist, (p0x, p0y), 0.0
+
+            # project (p - p0) onto seg: t = dot(p-p0, seg) / |seg|^2
+            apx = px - p0x
+            apy = py - p0y
+            t = (apx * sx + apy * sy) / seg_len_sq
+
+            # clamp to segment
+            if t <= 0.0:
+                t_clamped = 0.0
+                cx, cy = p0x, p0y
+            elif t >= 1.0:
+                t_clamped = 1.0
+                cx, cy = p1x, p1y
+            else:
+                t_clamped = t
+                cx = p0x + sx * t_clamped
+                cy = p0y + sy * t_clamped
+
+            dx = px - cx
+            dy = py - cy
+            dist = math.hypot(dx, dy)
+            if dist <= snap_threshold:
+                distances[dist] = e
+
+        if distances:
+            min_dist = float("inf")
+            for key in distances:
+                if key < min_dist:
+                    min_dist = key
+
+            idx = snap_obj.obj.data.edges[distances[min_dist]].vertices
+            edge_verts = (snap_obj.verts_3d[idx[0]], snap_obj.verts_3d[idx[1]])
+            division_points = divide_vector(
+                edge_verts[0], edge_verts[1], 2
+            )  # TODO Make it work for different divisions
+            for division_point in division_points:
+                intersection = tool.Cad.point_on_edge(division_point, (ray_target, loc))
+                distance = (division_point - intersection).length
+                if distance < snap_threshold:
+                    snap_point = {
+                        "object": snap_obj.obj,
+                        "type": "Edge Center",
+                        "point": division_point.copy(),
+                        "distance": distance,
+                    }
+                    points.append(snap_point)
+
+            intersection = tool.Cad.intersect_edges_v2((ray_target, loc), edge_verts)
+            if intersection[0]:
+                if tool.Cad.is_point_on_edge(intersection[1], edge_verts):
+                    distance = (intersection[1] - intersection[0]).length
+                    if distance < snap_threshold:
+                        snap_point = {
+                            "object": snap_obj.obj,
+                            "type": "Edge",
+                            "point": intersection[1].copy(),
+                            "edge_verts": edge_verts,
+                            "distance": distance,
+                        }
+                        points.append(snap_point)
+
+        return points
+
     @classmethod
     def ray_cast_by_proximity(
         cls,
@@ -457,7 +601,8 @@ class Raycast(bonsai.core.tool.Raycast):
             if bbox_2d:
                 if tool.Raycast.intersect_mouse_2d_bounding_box(mouse_pos, bbox_2d):
                     if tool.Raycast.object_is_visible_in_clipping_plane(obj):
-                        objs_to_raycast.append(obj)
+                        snap_obj = cls.create_snap_obj(obj)
+                        objs_to_raycast.append(snap_obj)
 
         return objs_to_raycast
 
@@ -474,12 +619,6 @@ class Raycast(bonsai.core.tool.Raycast):
         face_index = None
         # Wireframes
         if obj.type in {"EMPTY", "CURVE"} or (hasattr(obj.data, "polygons") and len(obj.data.polygons) == 0):
-            snap_points = tool.Raycast.ray_cast_by_proximity(context, event, obj)
-            if snap_points:
-                hit = sorted(snap_points, key=lambda x: x["distance"])[0]["point"]
-                if hit:
-                    hit_world = obj.original.matrix_world @ hit
-                    return obj, hit_world, face_index
             return None, None, None
         # Meshes
         else:
@@ -514,19 +653,20 @@ class Raycast(bonsai.core.tool.Raycast):
 
         ray_origin, ray_target, ray_direction = cls.get_viewport_ray_data(context, event)
 
-        for obj in objs_to_raycast:
+        for snap_obj in objs_to_raycast:
             if not include_wireframes and (
-                obj.type in {"EMPTY", "CURVE"} or (hasattr(obj.data, "polygons") and len(obj.data.polygons) == 0)
+                snap_obj.obj.type in {"EMPTY", "CURVE"}
+                or (hasattr(snap_obj.obj.data, "polygons") and len(snap_obj.obj.data.polygons) == 0)
             ):
                 continue
 
-            snap_obj, hit, face_index = cls.cast_rays_to_single_object(context, event, obj)
+            hit_obj, hit, face_index = cls.cast_rays_to_single_object(context, event, snap_obj.obj)
 
             if hit is not None:
                 length_squared = (hit - ray_origin).length_squared
                 if best_obj is None or length_squared < best_length_squared:
                     best_length_squared = length_squared
-                    best_obj = snap_obj
+                    best_obj = hit_obj
                     best_hit = hit
                     best_face_index = face_index
 
@@ -547,3 +687,261 @@ class Raycast(bonsai.core.tool.Raycast):
         if lens < 50:
             snap_threshold *= value
         return snap_threshold
+
+    @classmethod
+    def create_snap_obj(cls, obj):
+        # TODO Will this be done for every object or only wireframe?
+        for snap_obj in cls.snap_objs:
+            if obj.name == snap_obj.obj.name:
+                return snap_obj
+        snap_obj = SnapObj(obj)
+        cls.snap_objs.append(snap_obj)
+        return snap_obj
+
+    @classmethod
+    def clear_snap_objs(cls):
+        TreeNode.__clear_all__()
+        SnapObj.__clear_all__()
+        cls.snap_objs.clear()
+
+
+class TreeNode:
+    all = []
+
+    def __init__(self, box: tuple):
+        self.__class__.all.append(self)
+        self.box = box
+        self.child_a = None
+        self.child_b = None
+        self.edges = []
+
+    def __clear_all__():
+        for instance in TreeNode.all:
+            del instance
+        TreeNode.all.clear()
+
+
+class SnapObj:
+    max_depth = 9
+    all = []
+
+    def __init__(self, obj: bpy.types.Object):
+        self.__class__.all.append(self)
+        self.obj = obj
+        self.root = self._create_root_node()
+        self.root.edges = [e.index for e in obj.data.edges]
+        self.split_box(self.root, 0)
+        self.verts_3d = [obj.matrix_world @ v.co for v in obj.data.vertices]
+        self.snap_points = []
+
+    def __clear_all__():
+        for instance in SnapObj.all:
+            del instance
+        SnapObj.all.clear()
+
+    def _create_root_node(self) -> TreeNode:
+        bbox = tool.Blender.get_object_bounding_box(self.obj)
+        min_point = self.obj.matrix_world @ bbox["min_point"]
+        max_point = self.obj.matrix_world @ bbox["max_point"]
+        new_bbox = self.expand_bounding_box((min_point, max_point))
+        return TreeNode(new_bbox)
+
+    def divide_bounding_box_along_longest_axis(
+        self, min_pt: Vector, max_pt: Vector
+    ) -> Union[tuple[Vector, Vector], tuple[Vector, Vector]]:
+        """
+        Divide a bounding box into two equal parts along the axis with the longest dimension.
+
+        Args:
+            min_pt: The minimum point of the bounding box.
+            max_pt: The maximum point of the bounding box.
+
+        Returns:
+            list: A list of two tuples, each containing the minimum and maximum points of the divided boxes.
+        """
+
+        # Calculate the dimensions of the box
+        dx = max_pt.x - min_pt.x
+        dy = max_pt.y - min_pt.y
+        dz = max_pt.z - min_pt.z
+
+        # Determine the axis with the longest dimension
+        if dx >= dy and dx >= dz:
+            # Divide along the x-axis
+            mid_x = min_pt.x + dx / 2
+            box1 = (min_pt, Vector((mid_x, max_pt.y, max_pt.z)))
+            box2 = (Vector((mid_x, min_pt.y, min_pt.z)), max_pt)
+        elif dy >= dx and dy >= dz:
+            # Divide along the y-axis
+            mid_y = min_pt.y + dy / 2
+            box1 = (min_pt, Vector((max_pt.x, mid_y, max_pt.z)))
+            box2 = (Vector((min_pt.x, mid_y, min_pt.z)), max_pt)
+        else:
+            # Divide along the z-axis
+            mid_z = min_pt.z + dz / 2
+            box1 = (min_pt, Vector((max_pt.x, max_pt.y, mid_z)))
+            box2 = (Vector((min_pt.x, min_pt.y, mid_z)), max_pt)
+
+        return [box1, box2]
+
+    def expand_bounding_box(self, box: tuple[Vector, Vector], offset: float = 0.1) -> tuple[Vector, Vector]:
+        """
+        Expand a 3D bounding box by a given offset.
+
+        Args:
+            min_pt: The minimum point of the bounding box.
+            max_pt: The maximum point of the bounding box.
+            offset: The offset to expand the bounding box by.
+
+        Returns:
+            tuple: A tuple containing the new minimum and maximum points of the expanded bounding box.
+        """
+
+        min_pt, max_pt = box
+        # Calculate the new minimum and maximum points
+        new_min_pt = Vector((min_pt.x - offset, min_pt.y - offset, min_pt.z - offset))
+        new_max_pt = Vector((max_pt.x + offset, max_pt.y + offset, max_pt.z + offset))
+
+        return new_min_pt, new_max_pt
+
+    def split_box(self, parent: TreeNode, depth: int):
+        """
+        Splits the bounding box creating two child nodes to compose a BVH Tree recursively.
+
+        Args:
+            parent: the TreeNode instance that represents the parent node of a BVH Tree.
+            depth: the depth of the BVH Tree no be used in recursion.
+        """
+        if depth > self.max_depth:
+            return
+        box_a, box_b = self.divide_bounding_box_along_longest_axis(parent.box[0], parent.box[1])
+        parent.child_a = TreeNode(box_a)
+        parent.child_b = TreeNode(box_b)
+        edges_a = []
+        edges_b = []
+        for e in parent.edges:
+            verts_idx = [v for v in self.obj.data.edges[e].vertices]
+            verts_coords = []
+            for idx in verts_idx:
+                if idx < len(self.obj.data.vertices):
+                    verts_coords.append(self.obj.matrix_world @ self.obj.data.vertices[idx].co)
+            if self.line_intersects_box(verts_coords[0], verts_coords[1], parent.child_a.box):
+                edges_a.append(e)
+            if self.line_intersects_box(verts_coords[0], verts_coords[1], parent.child_b.box):
+                edges_b.append(e)
+        parent.child_a.edges = edges_a
+        parent.child_b.edges = edges_b
+        self.split_box(parent.child_a, depth + 1)
+        self.split_box(parent.child_b, depth + 1)
+
+    def raycast_box(
+        self, context: bpy.types.Context, event: bpy.types.Event, node: TreeNode, rays: tuple[Vector, Vector]
+    ) -> bool:
+        """
+        Raycast bounding box.
+
+        Args:
+            context: Blender context.
+            event: Blender event.
+            node: a TreeNode instance.
+            rays: tuple containing ray origin and ray direction
+
+        Returns:
+            True if hits the box or False otherwise.
+        """
+        box = node.box
+        min_v = box[0]
+        max_v = box[1]
+        t_min = 0.0
+        t_max = float("inf")
+        ray_origin, ray_dir = rays
+        inv_dir = Vector((1.0 / r if r != 0.0 else 1e32) for r in (ray_dir.x, ray_dir.y, ray_dir.z))
+        # X
+        tx1 = (min_v.x - ray_origin.x) * inv_dir[0]
+        tx2 = (max_v.x - ray_origin.x) * inv_dir[0]
+        tmin = min(tx1, tx2)
+        tmax = max(tx1, tx2)
+        # Y
+        ty1 = (min_v.y - ray_origin.y) * inv_dir[1]
+        ty2 = (max_v.y - ray_origin.y) * inv_dir[1]
+        tmin = max(tmin, min(ty1, ty2))
+        tmax = min(tmax, max(ty1, ty2))
+        # Z
+        tz1 = (min_v.z - ray_origin.z) * inv_dir[2]
+        tz2 = (max_v.z - ray_origin.z) * inv_dir[2]
+        tmin = max(tmin, min(tz1, tz2))
+        tmax = min(tmax, max(tz1, tz2))
+        return (tmax >= max(tmin, t_min)) and (tmin <= t_max)
+
+    def line_intersects_box(self, v1: mathutils.Vector, v2: mathutils.Vector, box: tuple) -> bool:
+        """
+        Check if a line segment intersects an axis-aligned bounding box (AABB).
+
+        Args:
+            v1: The first endpoint of the line segment as a mathutils.Vector.
+            v2: The second endpoint of the line segment as a mathutils.Vector.
+            box: A tuple containing the minimum and maximum points of the AABB, where each point is a mathutils.Vector.
+
+        Returns:
+            bool: True if the segment [v1, v2] intersects the AABB; otherwise, False.
+        """
+        bmin, bmax = box
+        dir = v2 - v1
+        tmin = 0.0
+        tmax = 1.0
+
+        for i in range(3):
+            if abs(dir[i]) < 1e-12:
+                # Line is parallel to slab. If origin not within slab -> no hit.
+                if v1[i] < bmin[i] or v1[i] > bmax[i]:
+                    return False
+            else:
+                ood = 1.0 / dir[i]
+                t1 = (bmin[i] - v1[i]) * ood
+                t2 = (bmax[i] - v1[i]) * ood
+                if t1 > t2:
+                    t1, t2 = t2, t1
+                if t1 > tmin:
+                    tmin = t1
+                if t2 < tmax:
+                    tmax = t2
+                if tmin > tmax:
+                    return False
+
+        # If any overlap in [0,1] exists, there's intersection
+        return (tmax >= 0.0) and (tmin <= 1.0)
+
+    def raycast_boxes(
+        self,
+        context: bpy.types.Context,
+        event: bpy.Types.Event,
+        node: TreeNode,
+        intersected: Union[TreeNode] = [],
+        rays: tuple[Vector, Vector] = (),
+    ) -> Union[TreeNode]:
+        """
+        Raycast bounding box subdivisions recursively.
+
+        Args:
+            context: Blender context.
+            event: Blender event.
+            node: a TreeNode instance.
+            intersected: list of intersected boxes to use in recursion.
+            rays: tuple containing ray origin and ray direction
+
+        Returns:
+            tuple: a list of TreeNode instances that represent the subdivided boxes hit by the ray cast.
+        """
+        if not node.child_a:
+            intersected.append(node)
+            return intersected
+
+        intersects_a = self.raycast_box(context, event, node.child_a, rays)
+        intersects_b = self.raycast_box(context, event, node.child_b, rays)
+        if intersects_a:
+            intersected = self.raycast_boxes(context, event, node.child_a, intersected, rays)
+
+        if intersects_b:
+            intersected = self.raycast_boxes(context, event, node.child_b, intersected, rays)
+
+        return intersected
diff --git a/src/bonsai/bonsai/tool/snap.py b/src/bonsai/bonsai/tool/snap.py
index 02f484d388..1c1181866d 100644
--- a/src/bonsai/bonsai/tool/snap.py
+++ b/src/bonsai/bonsai/tool/snap.py
@@ -391,8 +391,8 @@ class Snap(bonsai.core.tool.Snap):
         # Wireframes
         # For wireframe we have to get all the objects so we can further calculate edge intersection
         for snap_obj in objs_to_raycast:
-            if snap_obj.type in {"EMPTY", "CURVE"} or (snap_obj.type == "MESH" and len(snap_obj.data.polygons) == 0):
-                snap_points = tool.Raycast.ray_cast_by_proximity(context, event, snap_obj)
+            if snap_obj.obj.type in {"EMPTY", "CURVE"} or (snap_obj.obj.type == "MESH" and len(snap_obj.obj.data.polygons) == 0):
+                snap_points = tool.Raycast.ray_cast_by_proximity_2d(context, event, snap_obj)
                 if snap_points:
                     for point in snap_points:
                         point["group"] = "Wireframe"
@@ -403,7 +403,7 @@ class Snap(bonsai.core.tool.Snap):
         ):
             results = []
             for obj in objs_to_raycast:
-                results.append(tool.Raycast.cast_rays_to_single_object(context, event, obj))
+                results.append(tool.Raycast.cast_rays_to_single_object(context, event, snap_obj.obj))
         else:
             results = []
             results.append(tool.Raycast.cast_rays_and_get_best_object(context, event, objs_to_raycast))