#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
GraphVisualization
==================
Module providing tools to create meshes and visualize graphs.
"""
__author__ = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__ = 'GPLv3 - GNU General Public License v3 (see LICENSE)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__ = 'https://idisco.info'
__download__ = 'https://www.github.com/ChristophKirst/ClearMap2'
import numpy as np
import functools as ft
import vispy.util.transforms as trf
import ClearMap.ParallelProcessing.SharedMemoryManager as smm
import ClearMap.ParallelProcessing.ParallelTraceback as ptb
import ClearMap.Analysis.Curves.Resampling as crs
###############################################################################
# ## Mesh generation
###############################################################################
[docs]
def mesh(graph=None, method='tubes', **kwargs):
if method == 'tubes':
return mesh_tube(graph, **kwargs)
else:
ValueError(f'Method n{method} not valid!')
###############################################################################
# ## Graph mesh using tubes
###############################################################################
[docs]
def mesh_tube(graph=None,
coordinates=None, radii=None, indices=None,
vertex_colors=None, edge_colors=None,
n_tube_points=None, default_radius=1,
color=None, default_color=(0.8, 0.1, 0.1, 1.0),
processes=None, verbose=False):
"""Construct mesh from edge geometry of a graph."""
if graph is not None:
if graph.has_edge_geometry(): # FIXME: split into two functions
if isinstance(coordinates, str):
coordinates, indices = graph.edge_geometry(name=coordinates, return_indices=True, as_list=False)
else:
raise ValueError(f'Expected coordinates to by None or str, found {coordinates}!')
try:
name = radii if radii is not None else 'radii'
radii = graph.edge_geometry(name=name, return_indices=False, as_list=False)
except KeyError:
print(f'No radii found in the graph, using uniform radii = {default_radius}!')
radii = default_radius * np.ones(coordinates.shape[0]) # FIXME: dtype=float
else:
if isinstance(coordinates, str) and coordinates != 'coordinates':
coordinates = graph.vertex_property(name=coordinates)
else:
coordinates = graph.vertex_coordinates()
indices = graph.edge_connectivity().flatten()
coordinates = np.vstack(coordinates[indices])
try:
radii = graph.vertex_radii()
radii = radii[indices]
except KeyError:
print(f'No radii found in the graph, using uniform radii = {default_radius}!')
radii = default_radius * np.ones(coordinates.shape[0]) # FIXME: dtype=float
n_edges = graph.n_edges
indices = np.array([2 * np.arange(0, n_edges), 2 * np.arange(1, n_edges+1)]).T
if vertex_colors is not None:
connectivity = graph.edge_connectivity()
edge_colors = (vertex_colors[connectivity[:, 0]] + vertex_colors[connectivity[:, 1]]) / 2.0
vertices, faces, colors = mesh_tube_from_coordinates_and_radii(coordinates, radii, indices,
n_tube_points=n_tube_points,
edge_colors=edge_colors,
processes=None)
return vertices, faces, colors
[docs]
def mesh_tube_from_coordinates_and_radii(coordinates, radii, indices, n_tube_points=8, edge_colors=None,
processes=None, verbose=False):
"""Construct a mesh from the edge geometry of a graph."""
coordinates_hdl = smm.insert(coordinates)
radii_hdl = smm.insert(radii)
indices_hdl = smm.insert(indices)
func = ft.partial(_parallel_mesh, coordinates_hdl=coordinates_hdl, radii_hdl=radii_hdl, indices_hdl=indices_hdl, n_tube_points=n_tube_points, verbose=verbose)
argdata = np.arange(len(indices))
# process in parallel
pool = smm.mp.Pool(processes=processes)
results = pool.map(func, argdata)
pool.close()
pool.join()
smm.free(coordinates_hdl)
smm.free(radii_hdl)
smm.free(indices_hdl)
results = list(results)
n_results = len(results)
vertices = [np.reshape(r[0], (-1, 3)) for r in results]
n_indices = [len(v) for v in vertices]
n_indices = np.cumsum(n_indices)
n_indices = np.hstack([[0], n_indices])
if edge_colors is not None:
# colors = [len(v) * [c] for v, c in zip(vertices, edge_colors)]
edge_colors = np.array(edge_colors)
if edge_colors.size in (3, 4):
edge_colors = [edge_colors] * len(vertices)
colors = [np.tile(c, [len(v), 1]) for v, c in zip(vertices, edge_colors)] # TODO: check we always want this broadcasting
colors = np.concatenate(colors)
else:
colors = None
vertices = np.concatenate(vertices)
faces = np.concatenate([results[i][1] + n_indices[i] for i in range(n_results)])
return vertices, faces, colors
def _mesh(coordinates, radii, n_tube_points=15, dtype='uint32'):
n_coordinates = len(coordinates)
tangents, normals, binormals = _frenet_frames(coordinates)
# circular tube
v = np.arange(n_tube_points, dtype=float) / n_tube_points * 2 * np.pi
c = np.cos(v)
s = np.sin(v)
r = radii[:, np.newaxis]
n = normals * r
b = binormals * r
# grid shape (npoints, ntube, 3)
grid = coordinates[:, np.newaxis, :] + c[np.newaxis, :, np.newaxis] * n[:, np.newaxis, :] + s[np.newaxis, :, np.newaxis] * b[:, np.newaxis, :]
# construct the mesh
n_segments = n_coordinates - 1
jp = np.ones(n_segments*n_tube_points, dtype=int)
jp[n_tube_points * np.arange(n_segments, dtype=int) - 1] -= n_tube_points
i1 = np.arange(n_segments*n_tube_points, dtype=int)
i2 = i1 + n_tube_points
i3 = i2 + jp
i4 = i1 + jp
indices = np.array(np.vstack([np.array([i1, i2, i4]).T, np.array([i2, i3, i4]).T]), dtype=dtype)
return grid, indices
def _frenet_frames(coordinates):
"""Calculates and returns the tangents, normals and binormals for a chain of coordinates."""
n_points = len(coordinates)
epsilon = 0.0001
# compute tangent vectors for each segment
tangents = np.roll(coordinates, -1, axis=0) - np.roll(coordinates, 1, axis=0)
tangents[0] = coordinates[1] - coordinates[0]
tangents[-1] = coordinates[-1] - coordinates[-2]
tangents = (tangents.T / np.linalg.norm(tangents, axis=1)).T
# get initial normal and binormal
t = np.abs(tangents[0])
smallest = np.argmin(t)
normal = np.zeros(3, dtype=float)
normal[smallest] = 1
vec = np.cross(tangents[0], normal)
normals = np.zeros((n_points, 3))
normals[0] = np.cross(tangents[0], vec)
# compute change along trajectory
theta = np.arccos(np.clip(np.sum(tangents[:-1] * tangents[1:], axis=1), -1, 1))
vec = np.cross(tangents[:-1], tangents[1:])
nrm = np.linalg.norm(vec, axis=1)
# compute normal and binormal vectors along the path
for i in range(n_points-1):
normals[i+1] = normals[i]
if nrm[i] > epsilon:
v = vec[i] / nrm[i]
normals[i+1] = trf.rotate(-np.degrees(theta[i]), v)[:3, :3].dot(normals[i+1])
binormals = np.cross(tangents, normals)
return tangents, normals, binormals
@ptb.parallel_traceback
def _parallel_mesh(i, coordinates_hdl, radii_hdl, indices_hdl, n_tube_points=15, verbose=False):
coordinates = smm.get(coordinates_hdl)
radii = smm.get(radii_hdl)
start, end = smm.get(indices_hdl)[i]
coordinates = coordinates[start:end]
radii = radii[start:end]
if verbose:
if i % 1000 == 0:
print(f'Mesh calculation {i} / {len(smm.get(indices_hdl))}.')
return _mesh(coordinates, radii, n_tube_points)
###############################################################################
# ## Render graph using intrepolation
###############################################################################
[docs]
def interpolate_edge_geometry(graph, smooth = 5, order = 2,
points_per_pixel = 0.5,
processes = None, verbose = False):
"""Smooth center lines and radii of the edge geometry."""
if not graph.has_edge_geometry('coordinates'):
raise ValueError('Graph has no edge geometry!')
coordinates, indices = graph.edge_geometry('coordinates', return_indices=True, as_list=False)
radii = graph.edge_geometry('radii', as_list=False)
# prepare result arrays
# indices_interp = np.array([_n_points_per_edge(i[1]-i[0], points_per_pixel=points_per_pixel) for i in indices]);
# indices_interp = np.cumsum(indices_interp);
# indices_interp = np.hstack([[0], indices_interp])
# indices_interp = np.array([indices_interp[:-1], indices_interp[1:]]).T;
# coordinates_interp = np.zeros((indices_interp[-1,1], coordinates.shape[1]), dtype=float);
# radii_interp = np.zeros(indices_interp[-1,1], dtype=float);
# process in parallel
coordinates_hdl = smm.insert(coordinates)
radii_hdl = smm.insert(radii)
indices_hdl = smm.insert(indices)
# coordinates_interp_hdl = smm.insert(coordinates_interp)
# radii_interp_hdl = smm.insert(radii_interp)
# indices_interp_hdl = smm.insert(indices_interp)
func = ft.partial(_parallel_interpolate,
coordinates_hdl=coordinates_hdl, radii_hdl=radii_hdl, indices_hdl=indices_hdl,
# coordinates_interp_hdl=coordinates_interp_hdl, radii_interp_hdl=radii_interp_hdl, indices_interp_hdl=indices_interp_hdl,
smooth=smooth, order=order, points_per_pixel=points_per_pixel, verbose=verbose)
argdata = np.arange(len(indices))
pool = smm.mp.Pool(processes=processes)
results = pool.map(func, argdata)
pool.close()
pool.join()
smm.free(coordinates_hdl)
smm.free(radii_hdl)
smm.free(indices_hdl)
results = list(results)
# smm.free(coordinates_interp_hdl)
# smm.free(radii_interp_hdl)
# smm.free(indices_interp_hdl)
indices_interp = np.array([len(r[1]) for r in results])
indices_interp = np.cumsum(indices_interp)
indices_interp = np.hstack([[0], indices_interp])
indices_interp = np.array([indices_interp[:-1], indices_interp[1:]]).T
coordinates_interp = np.vstack([r[0] for r in results])
radii_interp = np.hstack([r[1] for r in results])
return coordinates_interp, radii_interp, indices_interp
@ptb.parallel_traceback
def _parallel_interpolate(i, coordinates_hdl, radii_hdl, indices_hdl,
# coordinates_interp_hdl, radii_interp_hdl, indices_interp_hdl,
smooth=5, order=2, points_per_pixel=0.5, verbose=False):
coordinates = smm.get(coordinates_hdl)
radii = smm.get(radii_hdl)
start, end = smm.get(indices_hdl)[i]
# coordinates_interp = smm.get(coordinates_interp_hdl)
# radii_interp = smm.get(radii_interp_hdl)
# start_interp,end_interp = smm.get(indices_interp_hdl)[i]
coordinates = coordinates[start:end]
radii = radii[start:end]
if verbose:
if i % 1000 == 0:
print(f'Mesh interpolation {i} / {len(smm.get(indices_hdl))}.')
# coordinates_interp[start:end], radii_interp[start:end]=
n_points = _n_points_per_edge(end-start, points_per_pixel=points_per_pixel)
# n_points = end_interp - start_interp
# coordinates_interp[start_interp:end_interp], radii_interp[start_interp:end_interp] =
return _interpolate_edge(coordinates, radii, n_points=n_points,
smooth=smooth, order=order)
def _n_points_per_edge(n_pixel, points_per_pixel):
n_points = int(np.ceil(n_pixel * points_per_pixel))
n_points = max(2, n_points)
return n_points
def _interpolate_edge(coordinates, radii, n_points, smooth=5, order=2):
order = min(coordinates.shape[0]-1, order)
coordinates_interp = crs.resample(coordinates, n_points=n_points, smooth=smooth, order=order)
radii_interp = crs.resample(radii, n_points=n_points, smooth=smooth, order=order)
return coordinates_interp, radii_interp
###############################################################################
# ## Tests
###############################################################################
def _test():
import numpy as np
import ClearMap.Tests.Files as tf
import ClearMap.Analysis.graphs.graph_processing 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.has_edge_geometry()
g.vertex_coordinates()
s = g.skeleton()
assert np.all(s == skeleton)
gc = gp.clean_graph(g, verbose=True)
gr = gp.reduce_graph(gc, verbose=True)
coordinates, indices = gr.edge_geometry(return_indices=True, as_list=False)
gr.set_edge_geometry(name='radii', values=np.ones(len(coordinates)))
radii = gr.edge_geometry('radii', as_list=False)
import ClearMap.Analysis.graphs.GraphVisualization as gv
reload(gv)
grid, grid_indices = gv.mesh_from_edge_geometry(coordinates=coordinates, radii=radii, indices=indices)
import ClearMap.Visualization.GraphVisual as gvi
gvi.plot_mesh(grid, grid_indices)
