diff --git a/src/deep_neurographs/feature_extraction.py b/src/deep_neurographs/feature_extraction.py
index 33aba7b..72c33e1 100644
--- a/src/deep_neurographs/feature_extraction.py
+++ b/src/deep_neurographs/feature_extraction.py
@@ -8,19 +8,23 @@
 
 """
 
-import numpy as np
 from copy import deepcopy
-from deep_neurographs import utils, geometry_utils
 from random import sample
+
+import numpy as np
 from scipy.linalg import svd
 
+from deep_neurographs import geometry_utils, utils
+
 NUM_IMG_FEATURES = 0
 NUM_SKEL_FEATURES = 9
 NUM_PC_FEATURES = 0
 
 
 # -- Wrappers --
-def generate_mutable_features(neurograph, img=True, pointcloud=True, skel=True):
+def generate_mutable_features(
+    neurograph, img=True, pointcloud=True, skel=True
+):
     features = dict()
     if img:
         features["img"] = generate_img_features(neurograph)
@@ -74,15 +78,15 @@ def generate_mutable_skel_features(neurograph):
         radius_i, radius_j = get_radii(neurograph, edge)
 
         dot1, dot2, dot3 = get_directionals(neurograph, edge, 5)
-        ddot1, ddot2, ddot3 = get_directionals(neurograph, edge, 10)
-        features[edge] = np.concatenate((length, dot1, dot2, dot3), axis=None)
+        ddot1, ddot2, ddot3 = get_directionals(neurograph, edge, 5)
+        features[edge] = np.concatenate((length, radius_i, radius_j, dot1, dot2, dot3), axis=None)
     return features
 
 
 def compute_length(neurograph, edge, metric="l2"):
     i, j = tuple(edge)
     xyz_1, xyz_2 = neurograph.get_edge_attr("xyz", i, j)
-    return utils.dist(xyz_1, xyz_2, metric=metric)
+    return geometry_utils.dist(xyz_1, xyz_2, metric=metric)
 
 
 def get_directionals(neurograph, edge, window_size):
@@ -91,15 +95,19 @@ def get_directionals(neurograph, edge, window_size):
     mutable_xyz_i, mutable_xyz_j = neurograph.get_edge_attr("xyz", i, j)
     mutable_xyz = np.array([mutable_xyz_i, mutable_xyz_j])
     mutable_tangent = geometry_utils.compute_tangent(mutable_xyz)
-    context_tangent_1 = geometry_utils.compute_context_vec(neurograph, i, mutable_tangent, window_size=window_size)
-    context_tangent_2 = geometry_utils.compute_context_vec(neurograph, j, mutable_tangent, window_size=window_size)
-    
+    context_tangent_i = geometry_utils.compute_context_vec(
+        neurograph, i, mutable_tangent, window_size=window_size
+    )
+    context_tangent_j = geometry_utils.compute_context_vec(
+        neurograph, j, mutable_tangent, window_size=window_size
+    )
+
     # Compute features
-    inner_product_1 = abs(np.dot(mutable_tangent, context_tangent_1))
-    inner_product_2 = abs(np.dot(mutable_tangent, context_tangent_2))
-    inner_product_3 = np.dot(context_tangent_1, context_tangent_2)
+    inner_product_1 = abs(np.dot(mutable_tangent, context_tangent_i))
+    inner_product_2 = abs(np.dot(mutable_tangent, context_tangent_j))
+    inner_product_3 = np.dot(context_tangent_i, context_tangent_j)
     return inner_product_1, inner_product_2, inner_product_3
-    
+
 
 def get_radii(neurograph, edge):
     i, j = tuple(edge)
@@ -114,15 +122,13 @@ def build_feature_matrix(neurographs, features, blocks):
     X = None
     block_to_idxs = dict()
     idx_to_edge = dict()
-    
+
     # Feature extraction
     for block_id in blocks:
         # Get features
         idx_shift = 0 if X is None else X.shape[0]
         X_i, y_i, idx_to_edge_i = build_feature_submatrix(
-            neurographs[block_id],
-            features[block_id],
-            idx_shift,
+            neurographs[block_id], features[block_id], idx_shift
         )
 
         # Concatenate
@@ -175,7 +181,6 @@ def combine_feature_vecs(features):
     return vec
 
 
-
 """
 
 def generate_node_features(neurograph, img=True, pointcloud=True, skel=True):
@@ -215,4 +220,4 @@ def generate_immutable_skel_features(neurograph):
 def _generate_immutable_skel_features(neurograph, edge):
     mean_radius = np.mean(neurograph.edges[edge]["radius"], axis=0)
     return np.concatenate((mean_radius), axis=None)
-"""
\ No newline at end of file
+"""
diff --git a/src/deep_neurographs/geometry_utils.py b/src/deep_neurographs/geometry_utils.py
index 582c088..d5121b3 100644
--- a/src/deep_neurographs/geometry_utils.py
+++ b/src/deep_neurographs/geometry_utils.py
@@ -1,17 +1,27 @@
 import numpy as np
-from deep_neurographs import utils
+from scipy.interpolate import CubicSpline, UnivariateSpline
 from scipy.linalg import svd
+from deep_neurographs import utils
 
 
 # Context Tangent Vectors
-def compute_context_vec(neurograph, i, mutable_tangent, window_size=5, return_pts=False, vec_type="tangent"):
+def compute_context_vec(
+    neurograph,
+    i,
+    mutable_tangent,
+    window_size=5,
+    return_pts=False,
+    vec_type="tangent",
+):
     # Compute context vecs
     branches = get_branches(neurograph, i)
     context_vec_list = []
     xyz_list = []
     ref_xyz = neurograph.nodes[i]["xyz"]
     for branch in branches:
-        context_vec, xyz = _compute_context_vec(branch, ref_xyz, window_size, vec_type)
+        context_vec, xyz = _compute_context_vec(
+            branch, ref_xyz, window_size, vec_type
+        )
         context_vec_list.append(context_vec)
         xyz_list.append(xyz)
 
@@ -69,10 +79,63 @@ def compute_svd(xyz):
 
 def compute_tangent(xyz):
     if xyz.shape[0] == 2:
-        tangent = (xyz[1] - xyz[0]) / utils.dist(xyz[1], xyz[0])
+        tangent = (xyz[1] - xyz[0]) / dist(xyz[1], xyz[0])
     else:
+        xyz = smooth_branch(xyz)
         U, S, VT = compute_svd(xyz)
         tangent = VT[0]
     return tangent / np.linalg.norm(tangent)
 
 
+# Smoothing
+def smooth_branch(xyz):
+    t = np.arange(len(xyz[:, 0]) + 12)
+    s = len(t) / 10
+    cs_x = UnivariateSpline(t, extend_boundary(xyz[:, 0]), s=s, k=3)
+    cs_y = UnivariateSpline(t, extend_boundary(xyz[:, 1]), s=s, k=3)
+    cs_z = UnivariateSpline(t, extend_boundary(xyz[:, 2]), s=s, k=3)
+    smoothed_x = trim_boundary(cs_x(t))
+    smoothed_y = trim_boundary(cs_y(t))
+    smoothed_z = trim_boundary(cs_z(t))
+    smoothed = np.column_stack((smoothed_x, smoothed_y, smoothed_z))
+    return smoothed
+
+
+def extend_boundary(x, num_boundary_points=6):
+    extended_x = np.concatenate(
+        (
+            np.linspace(x[0], x[1], num_boundary_points, endpoint=False),
+            x,
+            np.linspace(x[-2], x[-1], num_boundary_points, endpoint=False),
+        )
+    )
+    return extended_x
+
+
+def trim_boundary(x, num_boundary_points=6):
+    return x[num_boundary_points:-num_boundary_points]
+
+
+# Miscellaneous
+def compare_edges(xyx_i, xyz_j, xyz_k):
+    dist_ij = dist(xyx_i, xyz_j)
+    dist_ik = dist(xyx_i, xyz_k)
+    return dist_ij < dist_ik
+    
+
+def dist(x, y, metric="l2"):
+    """
+    Computes distance between "x" and "y".
+
+    Parameters
+    ----------
+
+    Returns
+    -------
+    float
+
+    """
+    if metric == "l1":
+        return np.linalg.norm(np.subtract(x, y), ord=1)
+    else:
+        return np.linalg.norm(np.subtract(x, y), ord=2)
\ No newline at end of file
diff --git a/src/deep_neurographs/graph_utils.py b/src/deep_neurographs/graph_utils.py
index 9fe321c..c2d4443 100644
--- a/src/deep_neurographs/graph_utils.py
+++ b/src/deep_neurographs/graph_utils.py
@@ -11,12 +11,23 @@
 
 from copy import deepcopy as cp
 
+import os
 import networkx as nx
 import numpy as np
 
 from deep_neurographs import swc_utils, utils
 
 
+def init_dense_graphs(swc_dir):
+    dense_graphs = dict()
+    for f in utils.listdir(swc_dir, ext=".swc"):
+        raw_txt = swc_utils.read_swc(os.path.join(swc_dir, f))
+        swc_dict = swc_utils.parse(raw_txt)
+        graph_id = f.replace(".0.swc", "")
+        graph = swc_utils.file_to_graph(swc_dict, graph_id=graph_id, set_attrs=True)
+        dense_graphs[graph_id] = graph
+    return dense_graphs
+
 def get_irreducibles(graph):
     leafs = []
     junctions = []
diff --git a/src/deep_neurographs/intake.py b/src/deep_neurographs/intake.py
index 9363f61..a748505 100644
--- a/src/deep_neurographs/intake.py
+++ b/src/deep_neurographs/intake.py
@@ -15,7 +15,7 @@
 from torch_geometric.data import Data
 
 from deep_neurographs import neurograph as ng
-from deep_neurographs import s3_utils, swc_utils, utils
+from deep_neurographs import geometry_utils, s3_utils, swc_utils, utils
 
 
 # --- Build graph ---
@@ -25,11 +25,12 @@ def build_neurograph(
     bucket=None,
     access_key_id=None,
     secret_access_key=None,
-    generate_mutables=True,
     max_mutable_degree=5,
     max_mutable_dist=50.0,
     prune=True,
     prune_depth=16,
+    origin=None,
+    shape=None,
 ):
     """
     Builds a neurograph from a directory of swc files, where each swc
@@ -37,7 +38,7 @@ def build_neurograph(
     other.
 
     """
-    neurograph = ng.NeuroGraph()
+    neurograph = ng.NeuroGraph(origin=origin, shape=shape)
     if bucket is not None:
         neurograph = init_immutables_from_s3(
             neurograph,
@@ -48,7 +49,6 @@ def build_neurograph(
             secret_access_key=secret_access_key,
             prune=prune,
             prune_depth=prune_depth,
-            smooth=True,
         )
     else:
         neurograph = init_immutables_from_local(
@@ -57,11 +57,10 @@ def build_neurograph(
             anisotropy=anisotropy,
             prune=prune,
             prune_depth=prune_depth,
-            smooth=True,
         )
     neurograph.generate_mutables(
-        max_degree=max_mutable_degree, max_dist=max_mutable_dist
-    )
+            max_degree=max_mutable_degree, max_dist=max_mutable_dist
+        )
     return neurograph
 
 
@@ -74,7 +73,7 @@ def init_immutables_from_s3(
     secret_access_key=None,
     prune=True,
     prune_depth=16,
-    smooth=True,
+    smooth=False,
 ):
     """
     To do...
diff --git a/src/deep_neurographs/neurograph.py b/src/deep_neurographs/neurograph.py
index 8a3ae79..8e3403d 100644
--- a/src/deep_neurographs/neurograph.py
+++ b/src/deep_neurographs/neurograph.py
@@ -18,11 +18,10 @@
 from more_itertools import zip_broadcast
 from plotly.subplots import make_subplots
 from scipy.spatial import KDTree
+from tifffile import imwrite
 
 from deep_neurographs import graph_utils as gutils
-from deep_neurographs import swc_utils, utils
-
-from tifffile import imwrite
+from deep_neurographs import geometry_utils, swc_utils, utils
 
 COLORS = list(mcolors.TABLEAU_COLORS.keys())
 nCOLORS = len(COLORS)
@@ -37,7 +36,7 @@ class NeuroGraph(nx.Graph):
 
     """
 
-    def __init__(self, label_mask=None):
+    def __init__(self, label_mask=None, origin=None, shape=None):
         """
         Parameters
         ----------
@@ -56,8 +55,17 @@ def __init__(self, label_mask=None):
         self.mutable_edges = set()
         self.target_edges = set()
         self.xyz_to_edge = dict()
+        if origin is not None and shape is not None:
+            self.init_bbox(origin, shape)
+        else:
+            self.bbox = None
 
     # --- Add nodes or edges ---
+    def init_bbox(self, origin, shape):
+        self.bbox = dict()
+        self.bbox["min"] = list(origin)
+        self.bbox["max"] = [self.bbox["min"][i] + shape[i] for i in range(3)]
+
     def generate_immutables(
         self, swc_id, swc_dict, prune=True, prune_depth=16
     ):
@@ -119,7 +127,7 @@ def generate_immutables(
         for j in junctions:
             self.junctions.add(node_id[j])
 
-    def generate_mutables(self, max_degree=5, max_dist=100.0):
+    def generate_mutables(self, max_degree=5, max_dist=50.0):
         """
         Generates edges for the graph.
 
@@ -133,21 +141,23 @@ def generate_mutables(self, max_degree=5, max_dist=100.0):
         self._init_kdtree()
         for leaf in self.leafs:
             xyz_leaf = self.nodes[leaf]["xyz"]
+            if not self.is_contained(xyz_leaf):
+                continue
             mutables = self._get_mutables(leaf, xyz_leaf, max_degree, max_dist)
             for xyz in mutables:
+                if not self.is_contained(xyz):
+                    continue
                 # Extract info on mutable connection
                 (i, j) = self.xyz_to_edge[xyz]
                 attrs = self.get_edge_data(i, j)
 
                 # Get connecting node
-                if utils.dist(xyz, attrs["xyz"][0]) < 5:
+                if geometry_utils.dist(xyz, attrs["xyz"][0]) < 5:
                     node = i
                     xyz = self.nodes[node]["xyz"]
-                elif utils.dist(xyz, attrs["xyz"][-1]) < 5:
+                elif geometry_utils.dist(xyz, attrs["xyz"][-1]) < 5:
                     node = j
                     xyz = self.nodes[node]["xyz"]
-                    if node == leaf:
-                        stop
                 else:
                     idxs = np.where(np.all(attrs["xyz"] == xyz, axis=1))[0]
                     node = self.add_immutable_node((i, j), attrs, idxs[0])
@@ -155,6 +165,7 @@ def generate_mutables(self, max_degree=5, max_dist=100.0):
                 # Add edge
                 self.add_edge(leaf, node, xyz=np.array([xyz_leaf, xyz]))
                 self.mutable_edges.add(frozenset((leaf, node)))
+        print("# mutable edges:", len(self.mutable_edges))
 
     def _get_mutables(self, query_id, query_xyz, max_degree, max_dist):
         """
@@ -170,22 +181,21 @@ def _get_mutables(self, query_id, query_xyz, max_degree, max_dist):
         None.
 
         """
-        # Search for connections
         best_xyz = dict()
         best_dist = dict()
         query_swc_id = self.nodes[query_id]["swc_id"]
         for xyz in self._query_kdtree(query_xyz, max_dist):
-            xyz = tuple(xyz) #.astype(int))
+            xyz = tuple(xyz)
             edge = self.xyz_to_edge[xyz]
             swc_id = gutils.get_edge_attr(self, edge, "swc_id")
             if swc_id != query_swc_id:
-                d = utils.dist(xyz, query_xyz)
-                if edge not in best_dist.keys():
-                    best_xyz[edge] = xyz
-                    best_dist[edge] = d
-                elif d < best_dist[edge]:
-                    best_xyz[edge] = xyz
-                    best_dist[edge] = d
+                d = geometry_utils.dist(xyz, query_xyz)
+                if swc_id not in best_dist.keys():
+                    best_xyz[swc_id] = xyz
+                    best_dist[swc_id] = d
+                elif d < best_dist[swc_id]:
+                    best_xyz[swc_id] = xyz
+                    best_dist[swc_id] = d
         return self._get_best_edges(best_dist, best_xyz, max_degree)
 
     def _get_best_edges(self, dist, xyz, max_degree):
@@ -257,67 +267,93 @@ def _init_kdtree(self):
         """
         self.kdtree = KDTree(list(self.xyz_to_edge.keys()))
 
-    def _query_kdtree(self, query, max_dist):
+    def _query_kdtree(self, query, dist):
         """
         Parameters
         ----------
         query : int
             Node id.
+        dist : float
+            Distance from query that is searched.
 
         Returns
         -------
         generator[tuple]
             Generator that generates the xyz coordinates cooresponding to all
-            nodes within a distance of "max_dist" from "query".
+            nodes within a distance of "dist" from "query".
 
         """
-        idxs = self.kdtree.query_ball_point(query, max_dist)
+        idxs = self.kdtree.query_ball_point(query, dist, return_sorted=True)
         return self.kdtree.data[idxs]
-
-
-    def init_targets(self, target_labels, log_path, anisotropy=[1.0, 1.0, 1.0], shift=[0, 0, 0]):
-        labels_to_edge = dict()
-        mistakes_xyz = utils.read_mistake_coords(log_path, anisotropy=anisotropy, shift=shift)
-        mistakes_kdtree = KDTree(mistakes_xyz)
+    
+    def init_targets_via_path(self, target_neurograph, target_densegraph):
+        site_to_site = dict()
+        pair_to_edge = dict()
         for edge in self.mutable_edges:
-            # Get edge
+            # Get projection
             i, j = tuple(edge)
-            xyz_1 = utils.to_img(self.nodes[i]["xyz"], anisotropy, shift=shift)
-            xyz_2 = utils.to_img(self.nodes[j]["xyz"], anisotropy, shift=shift)
-
-            # Check coords in bounds
-            bounds_1 = [xyz_1[i] >= target_labels.shape[i] for i in range(3)]
-            bounds_2 = [xyz_2[i] >= target_labels.shape[i] for i in range(3)]
-            if any(bounds_1) or any(bounds_2):
+            xyz_i, d_i = target_neurograph.get_projection(self.nodes[i]["xyz"])
+            xyz_j, d_j = target_neurograph.get_projection(self.nodes[j]["xyz"])
+            if d_i > 8 or d_j > 8:
                 continue
+            
+            # Check criteria
+            edge_i = target_neurograph.xyz_to_edge[xyz_i]
+            edge_j = target_neurograph.xyz_to_edge[xyz_j]
+            inclusion_i = xyz_i in site_to_site.keys()
+            inclusion_j = xyz_j in site_to_site.keys()
+            if edge_i != edge_j:
+                True
+            elif not inclusion_i and not inclusion_j:
+                site_to_site, pair_to_edge = self.add_site(
+                    site_to_site, pair_to_edge, xyz_i, xyz_j, edge
+                )
+            else:
+                xyz_i = xyz_i if inclusion_i else xyz_j
+                xyz_k = site_to_site[xyz_i] if inclusion_i else site_to_site[xyz_j]
+                if geometry_utils.compare_edges(xyz_i, xyz_j, xyz_k):
+                    site_to_site, pair_to_edge = self.remove_site(
+                        site_to_site, pair_to_edge, xyz_i, xyz_k
+                    )
+                    site_to_site, pair_to_edge = self.add_site(
+                        site_to_site, pair_to_edge, xyz_i, xyz_j, edge
+                    )
 
-            # Check img
-            cond_1 = utils.check_img_path(target_labels, xyz_1, xyz_2)
-
-            # Check mistake log
-            d_1, _ = mistakes_kdtree.query(xyz_1, k=1)
-            d_2, _ = mistakes_kdtree.query(xyz_2, k=1)
-            cond_2 = d_1 < 10 and d_2 < 10
-            if cond_1 and cond_2:
-                key = frozenset([self.nodes[i]["swc_id"], self.nodes[j]["swc_id"]])
-                if key in labels_to_edge.keys():
-                    # Check whether to update target edge
-                    cur_dist = self.compute_length(labels_to_edge[key])
-                    new_dist = self.compute_length(edge)
-                    if cur_dist < new_dist:
-                        continue
-                    else:
-                        self.target_edges.remove(labels_to_edge[key])
-
-                # Add new target edge
-                labels_to_edge[key] = edge
-                self.target_edges.add(edge)
         print("# target edges:", len(self.target_edges))
         print(
             "% target edges in mutable:",
             len(self.target_edges) / len(self.mutable_edges),
         )
 
+    def remove_site(self, site_to_site, pair_to_edge, xyz_i, xyz_j):
+        del site_to_site[xyz_i]
+        del site_to_site[xyz_j]
+        pair_to_edge = self._remove_pair_edge(pair_to_edge, xyz_i, xyz_j)
+        return site_to_site, pair_to_edge
+
+    def add_site(self, site_to_site, pair_to_edge, xyz_i, xyz_j, edge):
+        self.target_edges.add(edge)
+        site_to_site[xyz_i] = xyz_j
+        site_to_site[xyz_j] = xyz_i
+        pair_to_edge = self._add_pair_edge(pair_to_edge, xyz_i, xyz_j, edge)
+        return site_to_site, pair_to_edge
+
+    def _add_pair_edge(self, pair_to_edge, xyz_i, xyz_j, edge):
+        key = frozenset([xyz_i, xyz_j])
+        if key not in pair_to_edge.keys():
+            pair_to_edge[key] = set([edge])
+        else:
+            pair_to_edge[key].add(edge)
+        return pair_to_edge
+
+    def _remove_pair_edge(self, pair_to_edge, xyz_i, xyz_j):
+        key = frozenset([xyz_i, xyz_j])
+        edges = list(pair_to_edge[key])
+        if len(edges) == 1:
+            self.target_edges.remove(edges[0])
+            del pair_to_edge[key]
+        return pair_to_edge
+
     # --- Visualization ---
     def visualize_immutables(self, title="Immutable Graph", return_data=False):
         """
@@ -348,14 +384,19 @@ def visualize_mutables(self, title="Mutable Graph", return_data=False):
             return data
         else:
             utils.plot(data, title)
-        
 
-    def visualize_targets(self, target_graph=None, title="Target Edges", return_data=False):
+    def visualize_targets(
+        self, target_graph=None, title="Target Edges", return_data=False
+    ):
         data = [self._plot_nodes()]
         data.extend(self._plot_edges(self.immutable_edges, color="black"))
         data.extend(self._plot_edges(self.target_edges))
         if target_graph is not None:
-            data.extend(target_graph._plot_edges(target_graph.immutable_edges, color="blue"))
+            data.extend(
+                target_graph._plot_edges(
+                    target_graph.immutable_edges, color="blue"
+                )
+            )
         if return_data:
             return data
         else:
@@ -458,18 +499,33 @@ def num_mutables(self):
 
         """
         return len(self.mutable_edges)
-    
+
     def compute_length(self, edge, metric="l2"):
         i, j = tuple(edge)
         xyz_1, xyz_2 = self.get_edge_attr("xyz", i, j)
-        return utils.dist(xyz_1, xyz_2, metric=metric)
+        return geometry_utils.dist(xyz_1, xyz_2, metric=metric)
 
     def path_length(self, metric="l2"):
         length = 0
         for edge in self.immutable_edges:
             length += self.compute_length(edge, metric=metric)
         return length
-
+    
+    def get_projection(self, xyz):
+        _, idx = self.kdtree.query(xyz, k=1)
+        proj_xyz = tuple(self.kdtree.data[idx])
+        proj_dist = geometry_utils.dist(proj_xyz, xyz)
+        return proj_xyz, proj_dist
+
+    def is_contained(self, xyz):
+        if type(self.bbox) is dict:
+            for i in range(3):
+                lower_bool = xyz[i] < self.bbox["min"][i]
+                upper_bool = xyz[i] > self.bbox["max"][i]
+                if lower_bool or upper_bool:
+                    return False
+        return True
+                
     def get_edge_attr(self, key, i, j):
         attr_1 = self.nodes[i][key]
         attr_2 = self.nodes[j][key]
@@ -492,9 +548,63 @@ def to_line_graph(self):
         graph = nx.Graph()
         graph.add_nodes_from(self.nodes)
         graph.add_edges_from(self.edges)
-        return nx.line_graph(graph)        
+        return nx.line_graph(graph)
+        
+        
+"""
 
-    
+    def init_targets(
+        self,
+        target_labels,
+        log_path,
+        anisotropy=[1.0, 1.0, 1.0],
+        shift=[0, 0, 0],
+    ):
+        labels_to_edge = dict()
+        mistakes_xyz = utils.read_mistake_coords(
+            log_path, anisotropy=anisotropy, shift=shift
+        )
+        mistakes_kdtree = KDTree(mistakes_xyz)
+        for edge in self.mutable_edges:
+            # Get edge
+            i, j = tuple(edge)
+            xyz_1 = utils.to_img(self.nodes[i]["xyz"], anisotropy, shift=shift)
+            xyz_2 = utils.to_img(self.nodes[j]["xyz"], anisotropy, shift=shift)
+
+            # Check coords in bounds
+            bounds_1 = [xyz_1[i] >= target_labels.shape[i] for i in range(3)]
+            bounds_2 = [xyz_2[i] >= target_labels.shape[i] for i in range(3)]
+            if any(bounds_1) or any(bounds_2):
+                continue
+
+            # Check img
+            cond_1 = utils.check_img_path(target_labels, xyz_1, xyz_2)
+
+            # Check mistake log
+            d_1, _ = mistakes_kdtree.query(xyz_1, k=1)
+            d_2, _ = mistakes_kdtree.query(xyz_2, k=1)
+            cond_2 = d_1 < 10 and d_2 < 10
+            if cond_1 and cond_2:
+                key = frozenset(
+                    [self.nodes[i]["swc_id"], self.nodes[j]["swc_id"]]
+                )
+                if key in labels_to_edge.keys():
+                    # Check whether to update target edge
+                    cur_dist = self.compute_length(labels_to_edge[key])
+                    new_dist = self.compute_length(edge)
+                    if cur_dist < new_dist:
+                        continue
+                    else:
+                        self.target_edges.remove(labels_to_edge[key])
+
+                # Add new target edge
+                labels_to_edge[key] = edge
+                self.target_edges.add(edge)
+        print("# target edges:", len(self.target_edges))
+        print(
+            "% target edges in mutable:",
+            len(self.target_edges) / len(self.mutable_edges),
+        )
 
     def init_targets_old(self, log_path, dist_threshold):
         # Initializations
@@ -527,3 +637,5 @@ def init_targets_old(self, log_path, dist_threshold):
             "% target edges in mutable:",
             len(target_edges) / len(self.mutable_edges),
         )
+
+"""
\ No newline at end of file
diff --git a/src/deep_neurographs/swc_utils.py b/src/deep_neurographs/swc_utils.py
index 9738cf5..f2ab9d1 100644
--- a/src/deep_neurographs/swc_utils.py
+++ b/src/deep_neurographs/swc_utils.py
@@ -17,8 +17,7 @@
 import numpy as np
 from more_itertools import zip_broadcast
 
-from deep_neurographs import graph_utils as gutils
-from deep_neurographs import utils
+from deep_neurographs import geometry_utils, graph_utils as gutils, utils
 
 
 # -- io utils --
@@ -49,11 +48,7 @@ def parse(raw_swc, anisotropy=[1.0, 1.0, 1.0], idx=False):
 
     """
     # Initialize swc
-    swc_dict = {"id": [], "radius": [], "pid": []}
-    if idx:
-        swc_dict["idx"] = []
-    else:
-        swc_dict["xyz"] = []
+    swc_dict = {"id": [], "radius": [], "pid": [], "xyz": []}
 
     # Parse raw data
     min_id = np.inf
@@ -67,14 +62,9 @@ def parse(raw_swc, anisotropy=[1.0, 1.0, 1.0], idx=False):
             swc_dict["id"].append(int(parts[0]))
             swc_dict["radius"].append(float(parts[-2]))
             swc_dict["pid"].append(int(parts[-1]))
-            if idx:
-                swc_dict["idx"].append(
-                    read_idx(parts[2:5], anisotropy=anisotropy, offset=offset)
-                ) 
-            else:
-                swc_dict["xyz"].append(
-                    read_xyz(parts[2:5], anisotropy=anisotropy, offset=offset)
-                )
+            swc_dict["xyz"].append(
+                read_xyz(parts[2:5], anisotropy=anisotropy, offset=offset)
+            )
             if swc_dict["id"][-1] < min_id:
                 min_id = swc_dict["id"][-1]
 
@@ -109,29 +99,6 @@ def read_xyz(xyz, anisotropy=[1.0, 1.0, 1.0], offset=[0, 0, 0]):
     return tuple(xyz)
 
 
-def read_idx(xyz, anisotropy=[1.0, 1.0, 1.0], offset=[0, 0, 0]):
-    """
-    Reads the (z,y,x) coordinates from an swc file, then reverses and scales
-    them.
-
-    Parameters
-    ----------
-    zyx : str
-        (z,y,x) coordinates.
-    anisotropy : list[float]
-        Image to real-world coordinates scaling factors for [x, y, z] due to
-        anistropy of the microscope.
-
-    Returns
-    -------
-    list
-        The (x,y,z) coordinates from an swc file.
-
-    """
-    xyz = [int(float(xyz[i]) * anisotropy[i]) + offset[i] for i in range(3)]
-    return xyz
-
-
 def write_swc(path, contents):
     if type(content) is list:
         write_list(path, contents)
@@ -211,7 +178,7 @@ def write_graph(path, graph):
 
     """
     # loop through connected components
-    
+
     reindex = dict()
     edges = graph.edges if edge_list is None else edge_list
     for i, j in edges:
@@ -293,7 +260,7 @@ def smooth(swc_dict):
         graph = file_to_graph(swc_dict)
         leafs, junctions = gutils.get_irreducibles(graph)
         if len(junctions) == 0:
-            xyz = utils.smooth_branch(xyz)
+            xyz = geometry_utils.smooth_branch(xyz)
         else:
             idxs = []
             root = None
@@ -315,7 +282,7 @@ def smooth(swc_dict):
 
 def upd_edge(xyz, idxs):
     idxs = np.array(idxs)
-    xyz[idxs] = utils.smooth_branch(xyz[idxs].copy())
+    xyz[idxs] = geometry_utils.smooth_branch(xyz[idxs].copy())
     return xyz
 
 
diff --git a/src/deep_neurographs/utils.py b/src/deep_neurographs/utils.py
index de9c6c4..071e5e8 100644
--- a/src/deep_neurographs/utils.py
+++ b/src/deep_neurographs/utils.py
@@ -12,13 +12,11 @@
 import json
 import os
 import shutil
-import zarr
 
 import numpy as np
 import plotly.graph_objects as go
+import zarr
 from plotly.subplots import make_subplots
-from scipy.interpolate import UnivariateSpline, CubicSpline
-from scipy.linalg import svd
 
 
 # --- dictionary utils ---
@@ -114,6 +112,7 @@ def read_n5(path):
     """
     return zarr.open(zarr.N5FSStore(path), "r").volume
 
+
 def read_json(path):
     """
     Reads json file stored at "path".
@@ -239,24 +238,6 @@ def subplot(data1, data2, title):
 
 
 # --- miscellaneous ---
-def dist(x, y, metric="l2"):
-    """
-    Computes distance between "x" and "y".
-
-    Parameters
-    ----------
-
-    Returns
-    -------
-    float
-
-    """
-    if metric == "l1":
-        return np.linalg.norm(np.subtract(x, y), ord=1)
-    else:
-        return np.linalg.norm(np.subtract(x, y), ord=2)
-
-
 def pair_dist(pair_1, pair_2, metric="l2"):
     pair_1.reverse()
     d1 = _pair_dist(pair_1, pair_2)
@@ -272,35 +253,6 @@ def _pair_dist(pair_1, pair_2, metric="l2"):
     return max(d1, d2)
 
 
-def smooth_branch(xyz, round=True):
-    t = np.arange(len(xyz[:, 0]) + 12)
-    s = len(t) / 10
-    cs_x = UnivariateSpline(t, extend_boundary(xyz[:, 0]), s=s, k=3)
-    cs_y = UnivariateSpline(t, extend_boundary(xyz[:, 1]), s=s, k=3)
-    cs_z = UnivariateSpline(t, extend_boundary(xyz[:, 2]), s=s, k=3)
-    smoothed_x = trim_boundary(cs_x(t))
-    smoothed_y = trim_boundary(cs_y(t))
-    smoothed_z = trim_boundary(cs_z(t))
-    smoothed = np.column_stack((smoothed_x, smoothed_y, smoothed_z))
-    if round:
-        return smoothed #np.round(smoothed).astype(int)
-    else:
-        return smoothed
-
-
-def extend_boundary(x, num_boundary_points=6):
-    extended_x = np.concatenate((
-        np.linspace(x[0], x[1], num_boundary_points, endpoint=False),
-        x,
-        np.linspace(x[-2], x[-1], num_boundary_points, endpoint=False)
-    ))
-    return extended_x
-
-
-def trim_boundary(x, num_boundary_points=6):
-    return x[num_boundary_points:-num_boundary_points]
-
-
 def check_img_path(target_labels, xyz_1, xyz_2):
     d = dist(xyz_1, xyz_2)
     t_steps = np.arange(0, 1, 1 / d)
@@ -322,9 +274,8 @@ def check_img_path(target_labels, xyz_1, xyz_2):
                 return False
     ratio = len(collisions) / len(t_steps)
     return True if ratio > 1 / 3 else False
-        
-    
-    
+
+
 def line(xyz_1, xyz_2, t):
     return np.round((1 - t) * xyz_1 + t * xyz_2)