From 849c01678e3c90bb913e6a2c83830954e5669040 Mon Sep 17 00:00:00 2001 From: Forrest Collman Date: Wed, 18 Dec 2019 15:14:42 -0800 Subject: [PATCH 1/4] adding resample method --- meshparty/skeleton.py | 53 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 51 insertions(+), 2 deletions(-) diff --git a/meshparty/skeleton.py b/meshparty/skeleton.py index 0b13d55..f4c24f5 100644 --- a/meshparty/skeleton.py +++ b/meshparty/skeleton.py @@ -1,6 +1,6 @@ import numpy as np from meshparty import utils -from scipy import spatial, sparse +from scipy import spatial, sparse, interpolate from pykdtree.kdtree import KDTree as pyKDTree from copy import copy import json @@ -194,11 +194,60 @@ def cover_paths(self): @property def distance_to_root(self): - """ np.array : N length array with the distance to the root node along the skeleton. + """an an array of distances to root for each skeleton vertex + + Returns + ------- + np.array + N length array with the distance to the root node along the skeleton. """ ds = sparse.csgraph.dijkstra(self.csgraph, directed=False, indices=self.root) return ds + + def resample(self, spacing, kind='nearest'): + """ resample the skeleton to a new spacig + + Parameters + ---------- + spacing : [float] + desired edge spacing in units of vertices + """ + cpaths= self.cover_paths() + d_to_root = self.distance_to_root + path_counter=1 + branch_d = {self.root: 0} + vert_list= [np.array([self.vertices[self.root,:]])] + edge_list = [] + for path in cpaths: + cpath = path[:-1] + input_d = d_to_root[cpath] + des_d = np.arange(0,np.max(input_d), spacing) + fi = interpolate.interp1d(input_d, self.vertices[cpath,:], kind=kind) + new_verts = fi(des_d) + + # find the index of the old branch points in the new path + is_branch = np.isin(np.array(cpath), self.branch_points) + path_branch = cpath[is_branch] + path_branch_verts = self.vertices[path_branch, :] + tree = pyKDTree(new_verts) + new_branch_on_path = tree.query(path_branch_verts) + new_branch_on_path+=path_counter + + new_edges = np.vstack([np.arange(0, len(new_verts)-1), np.arange(1, len(new_verts))]).T + path_counter + new_branch_d = {pb: nw for pb, nw in zip(path_branch, new_branch_on_path)} + branch_d.update(new_branch_d) + end_node_new = branch_d[path[-1]] + new_edges = np.vstack([new_edges, [len(new_verts), end_node_new]]) + edge_list.append(new_edges) + vert_list.append(new_verts) + path_counter += len(new_verts) + + new_verts = np.vstack(vert_list) + new_edges = np.vstack(edge_list) + + return Skeleton(new_verts, new_edges, + root=0) def path_to_root(self, v_ind): ''' From c7c204d0e2cd425b4716a0fc9f690ea9e13e5f9d Mon Sep 17 00:00:00 2001 From: Forrest Collman Date: Thu, 19 Dec 2019 08:48:49 -0800 Subject: [PATCH 2/4] adding a test --- test/skeleton_test.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/skeleton_test.py b/test/skeleton_test.py index 92e2a63..1965f9e 100644 --- a/test/skeleton_test.py +++ b/test/skeleton_test.py @@ -109,4 +109,10 @@ def test_child_nodes(full_cell_skeleton): def test_downstream_nodes(full_cell_skeleton): sk = full_cell_skeleton assert len(sk.downstream_nodes(sk.root)) == sk.n_vertices - assert len(sk.downstream_nodes(300)) == 135 \ No newline at end of file + assert len(sk.downstream_nodes(300)) == 135 + +def test_resample(full_cell_skeleton): + sk = full_cell_skeleton + skn, resamp_map = sk.resample(1000) + assert(skn.n_vertices==4476) + assert(len(resamp_map)==sk.n_vertices) \ No newline at end of file From 9a21e82292f0a13a1717725f25c83d2f42e6fb0a Mon Sep 17 00:00:00 2001 From: Forrest Collman Date: Thu, 19 Dec 2019 08:48:56 -0800 Subject: [PATCH 3/4] resample function added --- meshparty/skeleton.py | 76 ++++++++++++++++++++++++++----------------- 1 file changed, 47 insertions(+), 29 deletions(-) diff --git a/meshparty/skeleton.py b/meshparty/skeleton.py index f4c24f5..911e9ff 100644 --- a/meshparty/skeleton.py +++ b/meshparty/skeleton.py @@ -206,48 +206,66 @@ def distance_to_root(self): return ds def resample(self, spacing, kind='nearest'): - """ resample the skeleton to a new spacig + """ resample the skeleton to a new spacing and filter it to only include components connected to root Parameters ---------- spacing : [float] desired edge spacing in units of vertices + + Returns + ------- + skeleton.Skeleton + a resampled skeleton, which has the parts which are not connected to root removed + resample_map: + a N long array with as many entries as vertices in the original skeleton. + the entry reflects which new skeleton vertex that vertex should be mapped to """ - cpaths= self.cover_paths() + cpaths= self.cover_paths d_to_root = self.distance_to_root - path_counter=1 - branch_d = {self.root: 0} - vert_list= [np.array([self.vertices[self.root,:]])] + path_counter=0 + branch_d = {} + vert_list= [] edge_list = [] + + resample_map = -np.ones(len(self.vertices), dtype=np.int32) + for path in cpaths: - cpath = path[:-1] - input_d = d_to_root[cpath] - des_d = np.arange(0,np.max(input_d), spacing) - fi = interpolate.interp1d(input_d, self.vertices[cpath,:], kind=kind) - new_verts = fi(des_d) - - # find the index of the old branch points in the new path - is_branch = np.isin(np.array(cpath), self.branch_points) - path_branch = cpath[is_branch] - path_branch_verts = self.vertices[path_branch, :] - tree = pyKDTree(new_verts) - new_branch_on_path = tree.query(path_branch_verts) - new_branch_on_path+=path_counter - - new_edges = np.vstack([np.arange(0, len(new_verts)-1), np.arange(1, len(new_verts))]).T + path_counter - new_branch_d = {pb: nw for pb, nw in zip(path_branch, new_branch_on_path)} - branch_d.update(new_branch_d) - end_node_new = branch_d[path[-1]] - new_edges = np.vstack([new_edges, [len(new_verts), end_node_new]]) - edge_list.append(new_edges) - vert_list.append(new_verts) - path_counter += len(new_verts) + if ~np.isinf(d_to_root[path[-1]]): + input_d = d_to_root[path] + des_d = np.arange(np.min(input_d),np.max(input_d), spacing) + fi = interpolate.interp1d(input_d, self.vertices[path,:], kind=kind, axis=0) + new_verts = fi(des_d) + + # find the index of the old branch points in the new path + is_branch = np.isin(np.array(path), self.branch_points) + path_branch = path[is_branch] + path_branch_verts = self.vertices[path_branch, :] + tree = pyKDTree(new_verts) + + map_ds, new_branch_on_path = tree.query(path_branch_verts) + map_ds, path_map = tree.query(self.vertices[path,:]) + resample_map[path]=path_map+path_counter + new_branch_on_path+=path_counter + + new_edges = np.vstack([ np.arange(len(new_verts)-1,0,-1), np.arange(len(new_verts)-2,-1,-1)]).T + path_counter + new_branch_d = {pb: nw for pb, nw in zip(path_branch, new_branch_on_path)} + branch_d.update(new_branch_d) + last_edge=self.edges[self.edges[:,0]==path[-1],:] + if len(last_edge)==1: + new_edges = np.vstack([new_edges, [path_counter, branch_d[last_edge[0,1]] ]]) + edge_list.append(new_edges) + vert_list.append(new_verts) + + path_counter += len(new_verts) + new_verts = np.vstack(vert_list) new_edges = np.vstack(edge_list) - + + return Skeleton(new_verts, new_edges, - root=0) + root=branch_d[self.root]), resample_map def path_to_root(self, v_ind): ''' From 50fa0f2632d1872d45a05ef0f3632ed3f0af30e6 Mon Sep 17 00:00:00 2001 From: Forrest Collman Date: Thu, 19 Dec 2019 08:57:55 -0800 Subject: [PATCH 4/4] adding comments --- meshparty/skeleton.py | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/meshparty/skeleton.py b/meshparty/skeleton.py index 911e9ff..8516b15 100644 --- a/meshparty/skeleton.py +++ b/meshparty/skeleton.py @@ -232,9 +232,13 @@ def resample(self, spacing, kind='nearest'): for path in cpaths: if ~np.isinf(d_to_root[path[-1]]): + # use the distance from root to parameterize the path input_d = d_to_root[path] + # the desired distances from root are evenly spaced according to spacing des_d = np.arange(np.min(input_d),np.max(input_d), spacing) + # setup an interpolation function based upon distance to root as input and xyz as output fi = interpolate.interp1d(input_d, self.vertices[path,:], kind=kind, axis=0) + # use the function to interpolate the new values new_verts = fi(des_d) # find the index of the old branch points in the new path @@ -242,28 +246,41 @@ def resample(self, spacing, kind='nearest'): path_branch = path[is_branch] path_branch_verts = self.vertices[path_branch, :] tree = pyKDTree(new_verts) - map_ds, new_branch_on_path = tree.query(path_branch_verts) - map_ds, path_map = tree.query(self.vertices[path,:]) - resample_map[path]=path_map+path_counter new_branch_on_path+=path_counter - - new_edges = np.vstack([ np.arange(len(new_verts)-1,0,-1), np.arange(len(new_verts)-2,-1,-1)]).T + path_counter + # create a temporary dictionary with this path's branch points new_branch_d = {pb: nw for pb, nw in zip(path_branch, new_branch_on_path)} + # update the overall mapping dictionary branch_d.update(new_branch_d) + + # map the entire path to the new vertices by euc distance + map_ds, path_map = tree.query(self.vertices[path,:]) + # update the mapping + resample_map[path]=path_map+path_counter + + # new edges just march down path from last vertex to first + new_edges = np.vstack([ np.arange(len(new_verts)-1,0,-1), np.arange(len(new_verts)-2,-1,-1)]).T + path_counter + # need to construct the last edge since it wasn't in the original + # find that last edge whose start point was the first vertex in the path last_edge=self.edges[self.edges[:,0]==path[-1],:] + # for the first path there won't be an edge as the first vertex is root if len(last_edge)==1: - new_edges = np.vstack([new_edges, [path_counter, branch_d[last_edge[0,1]] ]]) + # if we do have one, then we want to add an edge that is the from the first vertex + # in the path, to the new vertex (mapped through branch_d) of the edge we found + # in the original edge list + new_edges = np.vstack([new_edges, [path_counter, branch_d[last_edge[0,1]] ]]) + # collect the edges and vertices in a list edge_list.append(new_edges) vert_list.append(new_verts) - + # increment the counter to keep track of how many vertices we have path_counter += len(new_verts) - + # concatenate the results together new_verts = np.vstack(vert_list) new_edges = np.vstack(edge_list) - + # create a new skeleton + # TODO add options to update mesh mapping and vertex properties return Skeleton(new_verts, new_edges, root=branch_d[self.root]), resample_map