#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
GraphProcessing
===============
Module providing tools to create and manipulate graphs. Functions include graph
construction from skeletons, simplification, reduction and clean-up of graphs.
"""
__author__ = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__ = 'GPLv3 - GNU General Pulic License v3 (see LICENSE)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__ = 'http://idisco.info'
__download__ = 'http://www.github.com/ChristophKirst/ClearMap2'
import numpy as np
import ClearMap.IO.IO as io
import ClearMap.Analysis.Graphs.GraphGt as ggt;
import ClearMap.ImageProcessing.Topology.Topology3d as t3d
import ClearMap.ParallelProcessing.DataProcessing.ArrayProcessing as ap
import ClearMap.Utils.Timer as tmr;
###############################################################################
### Graphs from skeletons
###############################################################################
[docs]
def graph_from_skeleton(skeleton, points = None, radii = None, vertex_coordinates = True,
check_border = True, delete_border = False, verbose = False):
"""Converts a binary skeleton image to a graph-tool graph.
Arguments
---------
skeleton : array
Source with 2d/3d binary skeleton.
points : array
List of skeleton points as 1d indices of flat skeleton array (optional to save processing time).
radii : array
List of radii associated with each vertex.
vertex_coordinates : bool
If True, store coordiantes of the vertices / edges.
check_border : bool
If True, check if the boder is empty. The algorithm reuqires this.
delete_border : bool
If True, delete the border.
verbose : bool
If True, print progress information.
Returns
-------
graph : Graph class
The graph corresponding to the skeleton.
Note
----
Edges are detected between neighbouring foreground pixels using 26-connectivty.
"""
skeleton = io.as_source(skeleton);
if delete_border:
skeleton = t3d.delete_border(skeleton);
check_border = False;
if check_border:
if not t3d.check_border(skeleton):
raise ValueError('The skeleton array needs to have no points on the border!');
if verbose:
timer = tmr.Timer();
timer_all = tmr.Timer();
print('Graph from skeleton calculation initialize.!');
if points is None:
points = ap.where(skeleton.reshape(-1, order = 'A')).array;
if verbose:
timer.print_elapsed_time('Point list generation');
timer.reset();
#create graph
n_vertices = points.shape[0];
g = ggt.Graph(n_vertices=n_vertices, directed=False);
g.shape = skeleton.shape;
if verbose:
timer.print_elapsed_time('Graph initialized with %d vertices' % n_vertices);
timer.reset();
#detect edges
edges_all = np.zeros((0,2), dtype = int);
for i,o in enumerate(t3d.orientations()):
# calculate off set
offset = np.sum((np.hstack(np.where(o))-[1,1,1]) * skeleton.strides)
edges = ap.neighbours(points, offset);
if len(edges) > 0:
edges_all = np.vstack([edges_all, edges]);
if verbose:
timer.print_elapsed_time('%d edges with orientation %d/13 found' % (edges.shape[0], i+1));
timer.reset();
if edges_all.shape[0] > 0:
g.add_edge(edges_all);
if verbose:
timer.print_elapsed_time('Added %d edges to graph' % (edges_all.shape[0]));
timer.reset();
if vertex_coordinates:
vertex_coordinates = np.array(np.unravel_index(points, skeleton.shape, order=skeleton.order)).T;
g.set_vertex_coordinates(vertex_coordinates);
if radii is not None:
g.set_vertex_radius(radii);
if verbose:
timer_all.print_elapsed_time('Skeleton to Graph');
return g;
[docs]
def graph_to_skeleton(graph, sink=None, dtype=bool, values=True):
"""Create a binary skeleton from a graph."""
if not graph.has_edge_geometry():
coordinates = np.asarray(graph.vertex_coordinates(), dtype=int);
else:
coordinates = graph.edge_geometry('coordinates');
coordinates = np.asarray(np.vstack(coordinates), dtype=int);
if sink is None:
shape = graph.shape;
if shape is None:
shape = tuple(np.max(coordinates, axis = 0));
sink = np.zeros(shape, dtype=dtype);
coordinates = tuple(coordinates.T);
sink[coordinates] = values;
return sink;
###############################################################################
### Graph CleanUp
###############################################################################
[docs]
def mean_vertex_coordinates(coordinates):
return np.mean(coordinates, axis=0);
[docs]
def clean_graph(graph, remove_self_loops = True, remove_isolated_vertices = True,
vertex_mappings = {'coordinates' : mean_vertex_coordinates,
'radii' : np.max},
verbose = False):
"""Remove all cliques to get pure branch structure of a graph.
Arguments
---------
graph : Graph
The graph to clean up.
verbose : bool
If True, prin progress information.
Returns
-------
graph : Graph
A graph removed of all cliques.
Note
----
cliques are replaced by a single vertex connecting to all non-clique neighbours
The center coordinate is used for that vertex as the coordinate.
"""
if verbose:
timer = tmr.Timer();
timer_all = tmr.Timer();
# find branch points
n_vertices = graph.n_vertices;
degrees = graph.vertex_degrees();
branches = degrees >= 3;
n_branch_points = branches.sum();
if verbose:
timer.print_elapsed_time('Graph cleaning: found %d branch points among %d vertices' % (n_branch_points, n_vertices));
timer.reset();
if n_branch_points == 0:
return graph.copy();
# detect 'cliques', i.e. connected components of branch points
gb = graph.sub_graph(vertex_filter=branches, view=True);
components, counts = gb.label_components(return_vertex_counts=True);
#note: graph_tools components of a view is a property map of the full graph
components[branches] += 1;
# group components
components = _group_labels(components);
#note: graph_tools components of a view is a property map of the full graph
#note: remove the indices of the non branch nodes
components = components[1:];
clique_ids = np.where(counts > 1)[0];
n_cliques = len(clique_ids);
if verbose:
timer.print_elapsed_time('Graph cleaning: detected %d cliques of branch points' % n_cliques);
timer.reset();
# remove cliques
g = graph.copy();
g.add_vertex(n_cliques);
#mappings
properties = {};
mappings = {};
for k in vertex_mappings.keys():
if k in graph.vertex_properties:
mappings[k] = vertex_mappings[k];
properties[k] = graph.vertex_property(k);
vertex_filter = np.ones(n_vertices + n_cliques, dtype = bool)
for i,ci in enumerate(clique_ids):
vi = n_vertices + i;
cc = components[ci];
#remove clique vertices
vertex_filter[cc] = False;
# get neighbours
neighbours = np.hstack([graph.vertex_neighbours(c) for c in cc])
neighbours = np.setdiff1d(np.unique(neighbours), cc);
# connect to new node
g.add_edge([[n, vi] for n in neighbours]);
#map properties
for k in mappings.keys():
g.set_vertex_property(k, mappings[k](properties[k][cc]), vertex=vi);
if verbose and i+1 % 10000 == 0:
timer.print_elapsed_time('Graph cleaning: reducing %d / %d' % (i+1, n_cliques));
#generate new graph
g = g.sub_graph(vertex_filter=vertex_filter);
if remove_self_loops:
g.remove_self_loops();
if verbose:
timer.print_elapsed_time('Graph cleaning: removed %d cliques of branch points from %d to %d nodes and %d to %d edges' % (n_cliques, graph.n_vertices, g.n_vertices, graph.n_edges, g.n_edges));
timer.reset();
# remove isolated vertices
if remove_isolated_vertices:
non_isolated = g.vertex_degrees() > 0;
g = g.sub_graph(vertex_filter = non_isolated)
if verbose:
timer.print_elapsed_time('Graph cleaning: Removed %d isolated nodes' % np.logical_not(non_isolated).sum());
timer.reset();
del non_isolated;
if verbose:
timer_all.print_elapsed_time('Graph cleaning: cleaned graph has %d nodes and %d edges' % (g.n_vertices, g.n_edges));
return g;
def _group_labels(array):
"""Helper grouping labels in a 1d array."""
idx = np.argsort(array);
arr = array[idx];
dif = np.diff(arr);
res = np.split(idx, np.where(dif > 0)[0]+1);
return res
###############################################################################
### Grapch reduction
###############################################################################
[docs]
def reduce_graph(graph, vertex_to_edge_mappings = {'radii' : np.max},
edge_to_edge_mappings = {'length' : np.sum},
edge_geometry = True, edge_length = None,
edge_geometry_vertex_properties = ['coordinates', 'radii'],
edge_geometry_edge_properties = None,
return_maps = False, verbose = False):
"""Reduce graph by replacing all vertices with degree two."""
if verbose:
timer = tmr.Timer();
timer_all = tmr.Timer();
print('Graph reduction: initialized.');
#ensure conditions for precessing step are fullfiled
if graph.is_view:
raise ValueError('Cannot process on graph view, prune graph before graph reduction.');
if not np.all(np.diff(graph.vertex_indices())==1) or int(graph.vertex_iterator().next()) != 0:
raise ValueError('Graph vertices not ordered!')
#copy graph
g = graph.copy();
#find non branching points, i.e. vertices with deg 2
branch_points = g.vertex_degrees() == 2;
non_branch_ids = np.where(branch_points)[0];
branch_points = np.logical_not(branch_points)
branch_ids = np.where(branch_points)[0];
n_non_branch_points = len(non_branch_ids);
n_branch_points = len(branch_ids);
del branch_points;
if verbose:
timer.print_elapsed_time('Graph reduction: Found %d branching and %d non-branching nodes' % (n_branch_points, n_non_branch_points));
timer.reset();
#mappings
if edge_to_edge_mappings is None:
edge_to_edge_mappings = {};
if vertex_to_edge_mappings is None:
vertex_to_edge_mappings = {};
if edge_length:
if 'length' not in edge_to_edge_mappings.keys():
edge_to_edge_mappings['length'] = np.sum;
if 'length' not in g.edge_properties:
g.add_edge_property('length', np.ones(g.n_edges, dtype = float));
edge_to_edge = {};
edge_properties = {};
edge_lists = {};
for k in edge_to_edge_mappings.keys():
if k in g.edge_properties:
edge_to_edge[k] = edge_to_edge_mappings[k];
edge_properties[k] = g.edge_property_map(k);
edge_lists[k] = [];
edge_to_edge_keys = edge_to_edge.keys();
vertex_to_edge = {};
vertex_properties = {};
vertex_lists = {};
for k in vertex_to_edge_mappings.keys():
if k in g.vertex_properties:
vertex_to_edge[k] = vertex_to_edge_mappings[k];
vertex_properties[k] = g.vertex_property_map(k);
vertex_lists[k] = [];
vertex_to_edge_keys = vertex_to_edge.keys();
#edge geometry
calculate_vertex_to_vertex_map = False;
calculate_edge_to_edge_map = False;
calculate_edge_to_vertex_map = False;
if edge_geometry:
calculate_edge_to_vertex_map = True;
reduced_edge_to_vertex_map = [];
if return_maps:
calculate_vertex_to_vertex_map = True;
reduced_vertex_to_vertex_map = [];
calculate_edge_to_vertex_map = True;
reduced_edge_to_vertex_map = [];
calculate_edge_to_edge_map = True;
reduced_edge_to_edge_map = [];
if edge_geometry_edge_properties is not None:
calculate_edge_to_edge_map = True;
reduced_edge_to_edge_map = [];
if edge_geometry_vertex_properties is None:
edge_geometry_vertex_properties = [];
if edge_geometry_edge_properties is None:
edge_geometry_edge_properties = [];
#direct edges between branch points
edge_list = [];
checked = np.zeros(g.n_edges, dtype=bool);
for i,v in enumerate(branch_ids):
for e in g.vertex_edges_iterator(v):
if e.target().out_degree() != 2:
eid = g.edge_index(e);
if not checked[eid]:
checked[eid] = True;
vertex_ids = [int(e.source()), int(e.target())];
edge_list.append(vertex_ids);
for k in edge_to_edge_keys:
edge_lists[k].append(edge_to_edge[k]([edge_properties[k][e]]));
for k in vertex_to_edge_keys:
vv = vertex_properties[k];
vv = [vv[e.source()], vv[e.target()]];
vertex_lists[k].append(vertex_to_edge[k](vv));
if calculate_edge_to_vertex_map:
reduced_edge_to_vertex_map.append(vertex_ids);
if calculate_edge_to_edge_map:
reduced_edge_to_edge_map.append([e]);
if verbose and (i+1) % 250000 == 0:
timer.print_elapsed_time('Graph reduction: Scanned %d/%d branching nodes, found %d branches' % (i+1, n_branch_points, len(edge_list)));
if verbose:
timer.print_elapsed_time('Graph reduction: Scanned %d/%d branching nodes, found %d branches' % (i +1, n_branch_points, len(edge_list)));
#reduce non-direct edges
checked = np.zeros(g.n_vertices, dtype=bool);
for i,v in enumerate(non_branch_ids):
if not checked[v]:
# extract branch
checked[v] = True;
vertex = g.vertex(v);
vertices, edges, endpoints, isolated_loop = _graph_branch(vertex);
#print([int(vv) for vv in vertices], [[int(ee.source()), int(ee.target())] for ee in edges], endpoints, isolated_loop)
if not isolated_loop:
vertices_ids = [int(vv) for vv in vertices];
checked[vertices_ids[1:-1]] = True;
edge_list.append([int(ep) for ep in endpoints]);
#mappings
for k in edge_to_edge.keys():
ep = edge_properties[k];
edge_lists[k].append(edge_to_edge[k]([ep[e] for e in edges]))
for k in vertex_to_edge.keys():
vp = vertex_properties[k];
vertex_lists[k].append(vertex_to_edge[k]([vp[vv] for vv in vertices]));
#if edge_geometry is not None:
# branch_start.append(index_pos);
# index_pos += len(vertices_ids);
# branch_end.append(index_pos);
if calculate_edge_to_vertex_map:
reduced_edge_to_vertex_map.append(vertices_ids);
if calculate_edge_to_edge_map:
reduced_edge_to_edge_map.append(edges);
if verbose and (i+1) % 250000 == 0:
timer.print_elapsed_time('Graph reduction: Scanned %d/%d non-branching nodes found %d branches' % (i+1, n_non_branch_points, len(edge_list)));
if verbose:
timer.print_elapsed_time('Graph reduction: Scanned %d/%d non-branching nodes found %d branches' % (i+1, n_non_branch_points, len(edge_list)));
#reduce graph
#redefine branch edges
gr = g.sub_graph(edge_filter=np.zeros(g.n_edges, dtype=bool));
gr.add_edge(edge_list);
#determine edge ordering
edge_order = gr.edge_indices();
# remove non-branching points
if calculate_vertex_to_vertex_map:
gr.add_vertex_property('_vertex_id_', graph.vertex_indices());
gr = gr.sub_graph(vertex_filter=np.logical_not(checked));
if calculate_vertex_to_vertex_map:
reduced_vertex_to_vertex_map = gr.vertex_property('_vertex_id_');
gr.remove_vertex_property('_vertex_id_');
# maps
if calculate_edge_to_vertex_map:
reduced_edge_to_vertex_map = np.array(reduced_edge_to_vertex_map, dtype=object)[edge_order];
if calculate_edge_to_edge_map:
reduced_edge_to_edge_map = np.array(reduced_edge_to_edge_map, dtype=object)[edge_order];
# add edge properties
for k in edge_to_edge_keys:
gr.define_edge_property(k, np.array(edge_lists[k])[edge_order]);
for k in vertex_to_edge_keys:
gr.define_edge_property(k, np.array(vertex_lists[k])[edge_order]);
if edge_geometry is not None:
branch_indices = np.hstack(reduced_edge_to_vertex_map);
indices = [0] + [len(m) for m in reduced_edge_to_vertex_map];
indices = np.cumsum(indices);
indices = np.array([indices[:-1], indices[1:]]).T;
#branch_start = np.array(branch_start);
#branch_end = np.array(branch_end);
#vertex properties
#indices = np.array([branch_start, branch_end]).T;
#indices = indices[edge_order];
indices_use = indices;
for p in edge_geometry_vertex_properties:
if p in g.vertex_properties:
values = g.vertex_property(p);
values = values[branch_indices];
gr.set_edge_geometry(name=p, values=values, indices=indices_use);
indices_use = None;
for p in edge_geometry_edge_properties:
if p in g.edge_properties:
values = g.edge_property_map(p);
#there is one edge less than vertices in each reduced edge !
values = [[values[e] for e in edge] + [values[edge[-1]]] for edge in reduced_edge_to_edge_map];
gr.set_edge_geometry(name='edge_' + p, values=values, indices=indices_use);
indices_use = None;
if calculate_edge_to_edge_map:
reduced_edge_to_edge_map = [[g.edge_index(e) for e in edge] for edge in reduced_edge_to_edge_map]
if verbose:
timer_all.print_elapsed_time('Graph reduction: Graph reduced from %d to %d nodes and %d to %d edges' % (graph.n_vertices, gr.n_vertices, graph.n_edges, gr.n_edges));
if return_maps:
return gr, (reduced_vertex_to_vertex_map, reduced_edge_to_vertex_map, reduced_edge_to_edge_map);
else:
return gr;
def _graph_branch(v0):
"""Helper to follow a branch in a graph consisting of a chain of vertices of degree 2."""
edges_v0 = [e for e in v0.out_edges()];
assert(len(edges_v0) == 2)
vertices_1, edges_1, endpoint_1 = _follow_branch(edges_v0[0], v0);
if endpoint_1 != v0:
vertices_2, edges_2, endpoint_2 = _follow_branch(edges_v0[1], v0);
vertices_1.reverse();
vertices_1.append(v0);
vertices_1.extend(vertices_2);
edges_1.reverse();
edges_1.extend(edges_2);
isolated_loop = False;
else:
endpoint_2 = v0;
isolated_loop = True;
return (vertices_1, edges_1, [int(endpoint_1), int(endpoint_2)], isolated_loop);
def _follow_branch(e, v0):
"""Helper to follow branch from vertex v0 in direction of edge e."""
#argument v0 to prevent infinte loops for isolated closed loops
edges = [e];
v = e.target();
vertices = [v];
while v.out_degree() == 2 and v != v0:
for e_new in v.out_edges():
if e_new != e:
break;
e = e_new;
edges.append(e);
v = e.target();
vertices.append(v);
return (vertices, edges, v);
[docs]
def expand_graph_length(graph,length = 'length', return_edge_mapping = False):
"""Expand a reduced graph to a full graph."""
#data
ec = graph.edge_connectivity();
vertex_lengths = graph.edge_property(length);
edge_lengths = np.asarray(vertex_lengths, dtype=int) - 1;
#construct graph
new_edges = edge_lengths > 1; # edges to create
n_vertices = graph.n_vertices + int(np.sum(vertex_lengths - 2));
g = ggt.Graph(n_vertices=n_vertices);
#construct edges
offset = graph.n_vertices;
new_starts = np.hstack([[0], np.cumsum(edge_lengths[new_edges] - 1)]);
new_ends = new_starts[1:] - 1;
new_starts = new_starts[:-1];
new_to_new = np.array([np.arange(new_ends[-1]), np.arange(new_ends[-1])+1]).T
split = np.ones(new_ends[-1], dtype=bool).T;
split[new_ends[:-1]] = False;
new_to_new = new_to_new[split];
new_starts += offset;
new_ends += offset;
new_to_new += offset;
existing_to_existing = ec[np.logical_not(new_edges)];
existing_to_new = np.array([ec[new_edges,0], new_starts]).T;
new_to_existing = np.array([new_ends, ec[new_edges,1]]).T;
edges = np.vstack([existing_to_existing, existing_to_new, new_to_new, new_to_existing]);
g.add_edge(edges);
if return_edge_mapping:
#edge order in new graph
reorder = g.edge_indices();
# mapping[new_edge_id] = old_edge_id
existing_edges_id = np.where(np.logical_not(new_edges))[0]
new_edges_id = np.where(new_edges)[0];
mapping = np.hstack(
[existing_edges_id, #existing edges
new_edges_id, #existing to new
np.repeat(new_edges_id, new_ends - new_starts), #new to new
new_edges_id]); # new to existing;
mapping = mapping[reorder];
return g, mapping;
else:
return g;
[docs]
def expand_graph(graph):
raise NotImplementedError();
# if not graph.has_edge_geometry():
# return graph.copy();
#
# #calculate nodes of expanded graph
# lengths = graph.edge_geometry_lengths() - 2;
# n_add_vertices = np.sum(lengths);
#
# g = graph.copy();
# g.add_vertex(n_add_vertices);
# g.remove_edge_geometry();
#
# # re link edges and add vertex properties
# #current_vertex =
# #for e in graph.edge_iterator():
#
#
# #for name in graph.edge_geometry_names:
#
# graph.edge_geometry('coordinates');
###############################################################################
### Label Tracing
###############################################################################
[docs]
def trace_vertex_label(graph, vertex_label, condition, dilation_steps=1, max_iterations=None, pass_label=False, **condition_args):
"""Traces label within a graph.
Arguments
---------
vertex_label : array
Start label.
condition : function(graph, vertex)
A function determining if the vertex should be added to the labels.
steps : int
Number edges to jump to find new neighbours. Default is 1.
Returns
-------
vertex_label : array
The traced label.
"""
label = vertex_label.copy();
not_checked = np.logical_not(label);
if max_iterations is None:
max_iterations = np.inf;
iteration = 0;
while iteration < max_iterations:
label_dilated = graph.vertex_dilate_binary(label, steps=dilation_steps);
label_new = np.logical_xor(label, label_dilated);
label_new = np.logical_and(label_new, not_checked);
vertices_new = np.where(label_new)[0];
#print(label.sum(), label_dilated.sum(), label_new.sum())
#print(vertices_new)
if len(vertices_new) == 0:
break;
else:
if pass_label:
current_label = label.copy();
for v in vertices_new:
if condition(graph, v, current_label, **condition_args):
label[v] = True;
not_checked[v] = False;
else:
for v in vertices_new:
if condition(graph, v, **condition_args):
label[v] = True;
not_checked[v] = False;
iteration += 1;
return label;
[docs]
def trace_edge_label(graph, edge_label, condition, max_iterations = None, dilation_steps = 1, pass_label = False, **condition_args):
"""Traces label within a graph.
Arguments
---------
edge_label : array
Start label.
condition : function(graph, vertex)
A function determining if the vertex should be added to the labels.
steps : int
Number edges to jump to find new neighbours. Default is 1.
Returns
-------
edge_label : array
The traced label.
"""
edge_graph, edge_map = graph.edge_graph(return_edge_map=True)
traced = trace_vertex_label(edge_graph, edge_label, condition, max_iterations=max_iterations,
dilation_steps=dilation_steps, pass_label=pass_label, **condition_args)
label = np.zeros(len(traced), dtype=bool)
label[edge_map] = traced #TODO: check if initial label need to be reordered too ?
return label
# label = edge_label.copy();
# not_checked = np.logical_not(label);
# while True:
# label_dilated = graph.edge_dilate_binary(label, steps=steps);
# label_new = np.logical_xor(label, label_dilated);
# label_new = np.logical_and(label_new, not_checked);
# edges_new = np.where(label_new)[0];
# #print(label.sum(), label_dilated.sum(), label_new.sum())
# #print(vertices_new)
# if len(edges_new) == 0:
# break;
# else:
# for e in edges_new:
# if condition(graph, e):
# label[e] = True;
# not_checked[e] = False;
#
# return label;
###############################################################################
### Tests
###############################################################################
def _test():
import numpy as np
import ClearMap.Tests.Files as tf
import ClearMap.Analysis.Graphs.GraphProcessing as gp
#reload(gp)
skeleton = tf.source('skeleton');
#import ClearMap.Visualization.Plot3d as p3d
#p3d.plot(skeleton)
#reload(gp)
g = gp.graph_from_skeleton(skeleton)
g.vertex_coordinates()
s = gp.graph_to_skeleton(g)
assert np.all(s==skeleton)
gc = gp.clean_graph(g, verbose=True)
gr = gp.reduce_graph(gc, verbose=True)
gr2 = gr.copy();
gr2.set_edge_geometry_type('edge')
l = gr.edge_geometry_lengths();
print(l)
g = gp.ggt.Graph(n_vertices=10);
g.add_edge(np.array([[7,8],[7,9],[1,2],[2,3],[3,1],[1,4],[4,5],[2,6],[6,7]]));
g.set_vertex_coordinates(np.array([[10,10,10],[0,0,0],[1,1,1],[1,1,0],[5,0,0],[8,0,1],[0,7,1],[0,10,2],[0,12,3],[3,7,7]], dtype=float));
import ClearMap.Visualization.Plot3d as p3d
p3d.plot_graph_line(g)
gc = gp.clean_graph(g, verbose=True);
p3d.plot_graph_line(gc)
gr = gp.reduce_graph(gc, edge_geometry=True, verbose=True)
#gr.set_edge_geometry(0.1*np.ones(gr.edge_geometry(as_list=False).shape[0]), 'radii')
import ClearMap.Visualization.Plot3d as p3d
vertex_colors = np.random.rand(g.n_vertices, 4);
vertex_colors[:,3] = 1;
p3d.plot_graph_mesh(gr, default_radius=1, vertex_colors=vertex_colors)
eg = gr.edge_geometry(as_list=False);
egs = 0.5 * eg;
gr.set_edge_geometry(name='coordinates', values=egs)
#tracing
import numpy as np
import ClearMap.Visualization.Plot3d as p3d
import ClearMap.Analysis.Graphs.GraphProcessing as gp
g = gp.ggt.Graph(n_vertices=10);
g.add_edge(np.array([[0,1],[1,2],[2,3],[3,4],[4,0],[0,5],[5,6],[6,7],[0,8],[8,9],[9,0]]));
g.set_vertex_coordinates(np.array([[0,0,0],[1,0,0],[2,0,0],[2,2,0],[0,1,0],[0,0,1],[0,0,2],[0,0,3],[0,-1,0],[0,-1,-1]], dtype=float));
g.set_vertex_radii(np.array([10,6,4,6,7,4,2,2,5,5]) * 0.02)
vertex_colors = np.array([g.vertex_radii()]*4).T; vertex_colors = vertex_colors / vertex_colors.max(); vertex_colors[:,3] = 1;
p3d.plot_graph_mesh(g, default_radius=1, vertex_colors=vertex_colors)
def condition(graph, vertex):
r = graph.vertex_radii(vertex=vertex);
print('condition, vertex=%d, radius=%r' % (vertex,r))
return r >= 5 * 0.02;
label = np.zeros(g.n_vertices, dtype=bool);
label[0] = True;
traced = gp.trace_vertex_label(g, label, condition=condition, steps=1)
print(traced)
vertex_colors = np.array([[1,0,0,1],[1,1,1,1]])[np.asarray(traced, dtype=int)]
p3d.plot_graph_mesh(g, default_radius=1, vertex_colors=vertex_colors)
from importlib import reload;
reload(gp)
# edge tracing
import ClearMap.Analysis.Graphs.GraphGt as ggt
edges = [[0,1],[1,2],[2,3],[4,5],[5,6],[1,7]];
g = ggt.Graph(edges=edges)
l,m = g.edge_graph(return_edge_map=True);
import numpy as np
label = np.zeros(len(edges), dtype=bool);
label[1] = True;
import ClearMap.Analysis.Graphs.GraphProcessing as gp
def condition(graph, edge):
print('condition, edge=%d' % (edge,))
return True;
traced = gp.trace_edge_label(g, label, condition=condition);
print(traced)
# expansion of edge lengths
import numpy as np
import ClearMap.Tests.Files as tf
import ClearMap.Analysis.Graphs.GraphProcessing as gp
graph = gp.ggt.Graph(n_vertices=5);
graph.add_edge([[0,1],[0,2],[0,3],[2,3],[2,1],[0,4]]);
graph.add_edge_property('length', np.array([0,1,2,3,4,5])+2);
e, m = gp.expand_graph_length(graph, 'length', True)
import graph_tool.draw as gd
pos = gp.ggt.vertex_property_map_to_python(gd.sfdp_layout(e.base))
import matplotlib.pyplot as plt
plt.figure(1); plt.clf();
import matplotlib.collections as mc
import ClearMap.Visualization.Color as col
colors = np.array([col.color(c) for c in ['red', 'blue', 'green', 'black', 'purple', 'orange']])
ec = e.edge_connectivity();
lines = pos[ec];
cols = colors[m];
lc = mc.LineCollection(lines, linewidths=1, color=cols);
ax = plt.gca();
ax.add_collection(lc)
ax.autoscale()
ax.margins(0.1)
plt.scatter(pos[:,0], pos[:,1])
p =pos[:graph.n_vertices]
plt.scatter(p[:,0],p[:,1], color='red')
