import os
from pathlib import Path
from tqdm import tqdm
from ClearMap.Analysis.vasculature.geometry_utils import cartesian_to_polar, angle_between_vectors, compute_grid, \
interpolate_vectors
from ClearMap.Analysis.vasculature.vasc_graph_utils import set_artery_vein_if_missing, combine_arteries_and_veins, \
vertex_to_edge_property, edge_to_vertex_property, filter_graph_degrees, \
parallel_get_vessels_lengths, graph_gt_to_igraph, n_kinds, get_vertex_coordinates, vertex_filter_to_edge_filter
CORRECTED_GRAPH_BASE_NAME = 'data_graph_correcteduniverse'
try:
import cPickle as pickle
except ImportError: # python 3.x
import pickle
import numpy as np
import igraph
import ClearMap.IO.IO as clearmap_io
import ClearMap.Alignment.Annotation as annotation
import ClearMap.Analysis.Graphs.GraphGt as graph_gt
from ClearMap.Analysis.vasculature.flow.linear_system import LinearSystem
# TODO:
# - Change all orders to ids
# - Check that group_name is forced to controls whenever appropriate
# FIXME: check cases where group should be hard coded to controls (i.e. when it is the reference group)
N_PROCESSES = 15
VESSEL_PERMEABILITY_CONSTANT = 7.6 # TODO: check value
DEFAULT_UNITS = {'length': 'um', 'mass': 'ug', 'time': 'ms'}
[docs]
def load_graph(work_dir, group_name):
"""
Loads the graph from the given work directory and control.
Tries to load the default corrected graph, if it does not exist, tries to load the default graph.
"""
graph_dir = os.path.join(work_dir, group_name)
try:
graph = graph_gt.load(os.path.join(graph_dir, f'{CORRECTED_GRAPH_BASE_NAME}.gt'))
except FileNotFoundError:
graph = graph_gt.load(os.path.join(graph_dir, f'{group_name}_graph.gt'))
return graph
# FIXME: structures by name (here, cerebellum and hippocampus)
# FIXME: set to ID but won't work because algo seems designed for order
[docs]
def structure_specific_radius_filter(graph, artery_vein, label=None, specific_structures=((1006, 3), (463, 6))):
"""
This function reworks the graph so that the radius criterion is applied differently
depending on the structure (here, cerebellum and hippocampus)
Parameters
----------
graph
artery_vein
label
specific_structures
Returns
-------
"""
if label is None:
label = graph.vertex_annotation()
radii = graph.vertex_property('radii')
try:
arts = [] # @Sophie: what is arts ? arteries ? more explicit name
for id_, level in specific_structures:
label_leveled = annotation.convert_label(label, key='id', value='id', level=level)
structure = label_leveled == id_
arts.append(np.logical_and(radii <= 6, np.logical_and(structure, artery_vein))) # FIXME: magic numbers
artery_vein[structure] = 0
radii = graph.vertex_property('radii')
pb_art = arts[0]
for art in arts[1:]:
pb_art = np.logical_or(pb_art, art)
artery_vein = np.logical_or(pb_art, np.logical_and(artery_vein, radii >= 3)) # FIXME: magic numbers (6 in streamlines_flow.py)
except KeyError:
print(f'Could not find region, check json file')
return artery_vein
# WARNING: if clean_graph is True, This will use a structure specific radius filter
[docs]
def get_artery_vein_edge_coords_and_vectors(graph, mode, artery, vein, min_radius,
label=None, clean_graph=True,
distance_to_surface=None, min_dist_to_surface=2, invert=False,
coordinates_name='coordinates', orientation_criterion='distance',
return_filter=False):
artery_vein = combine_arteries_and_veins(graph, artery, vein, mode, min_radius)
if clean_graph:
artery_vein = structure_specific_radius_filter(graph, artery_vein, label=label)
if distance_to_surface is not None:
artery_vein = np.asarray(np.logical_and(artery_vein, distance_to_surface >= min_dist_to_surface))
vertex_filter = artery_vein
if invert: # i.e. capillaries
vertex_filter = np.logical_not(vertex_filter)
sub_graph = graph.sub_graph(vertex_filter=vertex_filter) # TODO: see why edge_connectivity is different for view (this would speedup otherwise)
out = graph_to_edge_coords_and_vectors(sub_graph, coordinates_name=coordinates_name,
orientation_criterion=orientation_criterion)
if return_filter:
out = list(out) + [vertex_filter_to_edge_filter(graph, vertex_filter)]
return out
[docs]
def compute_flow_f_schmid(work_dir, graph, cont, graph_name='correcteduniverse', dest_base_name='same'):
"""
Compute the flow and velocity of the graph using the method published in
Schmid F. et al. https://doi.org/10.1371/journal.pcbi.1005392
First, make an igraph version of the graph with the necassary properties
Then add:
- length (computed from the geometry of the edges)
- boundary_cap (True if edge is a boundary capillary)
- nkind (node kind, i.e. integer representing the type of node from
(0, 1, 2, 3) for (capillaries, arteries, veins, universe)
- pBC (pressure boundary conditions)
Parameters
----------
work_dir : str
The path to the work directory
graph: GraphGt.Graph
The graph to compute the flow and velocity from
cont : str
The name of the control group
graph_name : str
The name of the graph file
dest_base_name : str
The base name of the destination file. If 'same', the name will be the same as the graph name
Returns
-------
flow, velocity, pressure : ndarray, ndarray, ndarray
"""
tmp_graph = graph_gt_to_igraph(graph)
vessels_lengths = parallel_get_vessels_lengths(graph, n_processes=N_PROCESSES)
tmp_graph.es["length"] = vessels_lengths
graph_ml_path = os.path.join(work_dir, cont, f'{CORRECTED_GRAPH_BASE_NAME}.GraphML') # REFACTORING: could use path from graph ?
tmp_graph.write_graphml(graph_ml_path)
tmp_graph = igraph.read(graph_ml_path) # FIXME: is the read necessary
tmp_graph['defaultUnits'] = DEFAULT_UNITS
# Change capillary deadends adjacent to A or V to (n_kinds['arteries'], n_kinds['veins'])
tmp_graph.es['boundary_cap'] = [0] * tmp_graph.ecount()
for v in tmp_graph.vs(_degree_eq=1, nkind_eq=n_kinds['capillaries']):
tmp_graph.es[tmp_graph.incident(v)[0]]['boundary_cap'] = 1 # 0 by default, 1 for (degree1 and capillary)
# For boundary capillaries with diameter > 10, if vertices are connected to A or V, change to A or V
for e in tmp_graph.es(boundary_cap_eq=1, diameter_gt=10): # FIXME: magic number
for n_kind in (n_kinds['arteries'], n_kinds['veins']): # @ Sophie: is it correct to be sequential ?
if n_kind in tmp_graph.vs[e.tuple]['nkind']:
e['nkind'] = n_kind
for v in e.tuple:
tmp_graph.vs[v]['nkind'] = n_kind
# pressure fit. pBC = pressure boundary conditions
p = np.array([2.59476905e-05, -8.91218311e-03, 1.06888397e+00, 2.04594862e+01]) # FIXME: extract and explain where this comes from
tmp_graph.vs['pBC'] = [None] * tmp_graph.vcount()
# for all degree1 vertices, if veins, pBC = 10, else pBC = np.polyval(p, diameter)
tmp_graph.vs(_degree_eq=1, nkind_eq=n_kinds['veins'])['pBC'] = (
[10] * len(tmp_graph.vs(_degree_eq=1, nkind_eq=n_kinds['veins'])))
for v in tmp_graph.vs(_degree_eq=1, nkind_ne=n_kinds['veins']):
diameter = tmp_graph.es[tmp_graph.incident(v)[0]]['diameter']
v['pBC'] = np.polyval(p, diameter)
# p_a, _ = fit_for_pressureBC() # TODO: check if could remove
# assign_pressureBC(tmp_graph, p_a)
linear_system = LinearSystem(tmp_graph)
dest_base_name = graph_name if dest_base_name == 'same' else cont # REFACTOR: simplify logic
out_path = os.path.join(work_dir, cont, f'sampledict{dest_base_name}.pkl') # TODO: see if could extract path and share with load_sample_dict
sample_dict = linear_system.solve('iterative2', dest_file_path=out_path)
return [np.asarray(sample_dict[k]) for k in ['flow', 'v', 'pressure']]
[docs]
def compute_blood_flow(graph, work_dir, sample_name):
"""
Return the graph with the flow, velocity and pressure properties added
If this has not been computed before, it will be computed and saved to
the sample dict
Parameters
----------
graph : GraphGt.Graph
The graph to compute the flow, velocity and pressure from
work_dir : str
The path to the work directory
sample_name : str
The name of the sample (used in the file name)
Returns
-------
graph : GraphGt.Graph
The graph with the flow, velocity and pressure properties added
"""
flow, velocity, pressure = compute_flow_f_schmid(work_dir, graph, sample_name)
graph.add_edge_property('flow', flow)
graph.add_edge_property('veloc', velocity)
graph.add_vertex_property('pressure', pressure)
return graph
[docs]
def get_nb_radial_vessels(edge_color, radiality_threshold=0.7): # warning: unused
"""
Get the number of radial vessels from the edge color
Parameters
----------
edge_color : ndarray
radiality_threshold : float
Threshold to filter out vessels that are not sufficiently radial to the plane.
The default value is 0.7, which means that the vessel must be at least 70% radial to the plane to be
counted as radial. This is a good empirical value, but it may need to be adjusted for different
applications.
Returns
-------
nb_radial : int
The number of radial vessels
"""
radial = edge_color[:, 2] / edge_color.sum(axis=1) # TODO: check that edge_color.shape = (n, 3) and document
return np.sum(radial > radiality_threshold)
[docs]
def get_nb_parallel_vessels(edge_color, verbose=True, parallel_threshold=0.7): # warning: unused
"""
Get the number of parallel vessels from the edge color
Parameters
----------
edge_color : ndarray
parallel_threshold : float
Threshold to filter out vessels that are not sufficiently parallel to the plane.
The default value is 0.7, which means that the vessel must be at least 70% parallel to the plane to be
counted as parallel. This is a good empirical value, but it may need to be adjusted for different
applications.
verbose : bool
Whether to print the shape of the array
Returns
-------
nb_parallel : int
The number of parallel vessels
"""
planar = (edge_color[:, 0] + edge_color[:, 1]) / edge_color.sum(
axis=1) # TODO: check that edge_color.shape = (n, 3) and document
if verbose:
print(planar.shape)
return np.sum(planar > parallel_threshold)
[docs]
def graph_to_edge_coords_and_vectors(graph, coordinates_name='coordinates', orientation_criterion='distance'):
"""
Get the coordinates and vectors of the edges of the graph.
The edges can be oriented based on the criterion so that the vectors will follow this orientation.
Parameters
----------
graph : GraphGt.Graph
The graph to get the edge coordinates and vectors from
coordinates_name : str
The coordinates space to use. Typically, 'coordinates' or 'coordinates_atlas'
IT should be an existing vertex property of the graph
(that represents the coordinates of the vertices in some space)
orientation_criterion : str
The criterion to use to orient the edges. Can be either 'pressure', 'distance_to_surface' or '' (no sorting)
Returns
-------
"""
vertex_coordinates = get_vertex_coordinates(graph, coordinates_name)
connectivity = sort_connectivity_based_on_vertex_property(graph, orientation_criterion)
# Compute the edge coordinates as the mean of the coordinates of the vertices of the edge
edge_coordinates = np.round((vertex_coordinates[connectivity[:, 0]] +
vertex_coordinates[connectivity[:, 1]])
/ 2)
edge_vectors = coordinates_to_edge_vectors(vertex_coordinates, connectivity, normalized=False)
return edge_coordinates, edge_vectors
[docs]
def coordinates_to_edge_vectors(vertex_coordinates, connectivity, normalized=False):
"""
Get the edge vectors from the coordinates and connectivity.
The vectors are computed as the difference between the coordinates of the vertices of the edge.
.. warning::
If normalized is True, the vectors will be normalized to unit length.
If a normalized vector is 0, it will be replaced by a vector of 0s.
Parameters
----------
vertex_coordinates : ndarray
The coordinates of the vertices
connectivity : ndarray
The connectivity of the graph.
normalized : bool
Whether to normalize the vectors to unit length or not
Returns
-------
edge_vectors : ndarray
The vectors of the edges
"""
edges_coordinates = vertex_coordinates[connectivity]
edge_vectors = edges_coordinates[:, 1, :] - edges_coordinates[:, 0, :]
if normalized:
norms = np.linalg.norm(edge_vectors, axis=1)[:, None]
edge_vectors = np.divide(edge_vectors, norms, out=np.zeros_like(edge_vectors), where=norms != 0)
return edge_vectors
[docs]
def get_edge_vectors(graph, normed=False, orientation_criterion='distance',
coordinates_name='coordinates'):
"""
Get the edge vectors and sorted connectivity from the graph.
The edges can be oriented based on the criterion so that the vectors will follow this orientation.
If no criterion is given (''), the edges will not be oriented.
Parameters
----------
graph : GraphGt.Graph
The graph to get the edge vectors from
normed : bool
Whether to normalize the vectors to unit length or not
orientation_criterion : str
The criterion to use to orient the edges.
Can be either 'pressure', 'distance_to_surface' or '' (no sorting)
coordinates_name : str
The coordinates space to use. Typically, 'coordinates', 'coordinates_atlas',
'coordinates_scaled' or 'coordinates_MRI'
Returns
-------
"""
connectivity = sort_connectivity_based_on_vertex_property(graph, orientation_criterion)
vertex_coordinates = get_vertex_coordinates(graph, coordinates_name)
edge_vectors = coordinates_to_edge_vectors(vertex_coordinates, connectivity, normed)
return edge_vectors, connectivity
# FIXME: what are x, y, center and where do they come from ?
[docs]
def get_edge_vectors_spherical(graph, x, y, center, orientation_criterion='distance', normed=False):
connectivity = sort_connectivity_based_on_vertex_property(graph, orientation_criterion, reversed=True)
spherical_coord = np.array(cartesian_to_polar(x - center[0], y - center[1])).T
edge_vectors = coordinates_to_edge_vectors(spherical_coord, connectivity, normed)
return edge_vectors, connectivity
# FIXME: set criterion defaults but check usage first
[docs]
def sort_connectivity_based_on_vertex_property(graph, orientation_criterion, reversed=False):
"""
For each edge, sort the vertices based on the criterion.
i.e. get the edges direction based on the criterion,
which can be either 'pressure', 'distance_to_surface' or '' (no sorting)
Parameters
----------
graph
orientation_criterion : str
The criterion to use to sort the edges. Can be either 'pressure', 'distance_to_surface' or '' (no sorting)
Returns
-------
"""
if orientation_criterion.startswith('distance'):
orientation_criterion = 'distance_to_surface'
if orientation_criterion == '':
vertex_property = None
else:
try:
vertex_property = graph.vertex_property(orientation_criterion)
except KeyError as err:
print(f'Unknown criteria {orientation_criterion}, typical values are ("pressure", "distance")')
raise err
connectivity = graph.edge_connectivity()
if vertex_property is not None and vertex_property.shape[0] != 0:
sorted_indices = np.argsort(vertex_property[connectivity], axis=1)
if reversed:
sorted_indices = sorted_indices[:, ::-1]
# Apply the sorted indices to the connectivity
oriented_connectivity = connectivity[
np.arange(connectivity.shape[0])[:, None], sorted_indices] # TODO: TEST: compare to above
connectivity = oriented_connectivity
return connectivity
[docs]
def get_streamline_agv_path(work_dir, group_name, sample_name, mode):
return os.path.join(work_dir, f'streamline_AGV_{group_name}_{sample_name}_{mode}.npy')
[docs]
def get_streamline_aec_path(work_dir, group_name, sample_name, mode):
return os.path.join(work_dir, f'streamline_AEC_{group_name}_{sample_name}_{mode}.npy')
[docs]
def get_reference_streamlines(work_dir, reference_samples, mode):
AEC = []
AGV = []
for sample in reference_samples:
streamline_aec_path = get_streamline_aec_path(work_dir, 'controls', sample, mode)
streamline_agv_path = get_streamline_agv_path(work_dir, 'controls', sample, mode)
if not os.path.exists(streamline_aec_path) or not os.path.exists(streamline_agv_path):
take_avg_streamlines(work_dir, sample, mode=mode)
AEC.append(np.load(streamline_aec_path, allow_pickle=True))
AGV.append(np.load(streamline_agv_path, allow_pickle=True))
return AEC, AGV
[docs]
def take_avg_streamlines(work_dir, sample, mode='bigvessels', group_name='controls', min_radius=4,
coordinates_name='coordinates'):
graph = load_graph(work_dir, sample) # TODO: check if sample or group_name
graph = filter_graph_degrees(graph)
set_artery_vein_if_missing(graph)
artery = edge_to_vertex_property(graph, 'artery')
vein = edge_to_vertex_property(graph, 'vein')
d2s = graph.vertex_property('distance_to_surface')
art_edge_coordinates, art_graph_vector = (
get_artery_vein_edge_coords_and_vectors(graph, mode, artery, vein, min_radius=min_radius,
distance_to_surface=d2s, coordinates_name=coordinates_name))
np.save(get_streamline_aec_path(work_dir, group_name, sample, mode), art_edge_coordinates)
np.save(get_streamline_agv_path(work_dir, group_name, sample, mode), art_graph_vector)
[docs]
def avg_streamlines_grid(work_dir, reference_samples, mode='bigvessels', group_name='controls'): # FIXME: unused
AEC, AGV = get_reference_streamlines(work_dir, reference_samples, mode)
all_samples_flow_vectors = []
for arteries_edge_coordinates, arteries_graph_vector in zip(AEC, AGV):
grid = compute_grid(arteries_edge_coordinates)
flow_vectors = interpolate_vectors(arteries_graph_vector, arteries_edge_coordinates, grid)
all_samples_flow_vectors.append(flow_vectors)
flow_vectors_avgs = np.nanmean(np.array(all_samples_flow_vectors), axis=0)
dest_path = Path(work_dir) / f'streamline_grid_avg_{group_name}{mode}.npy'
np.save(dest_path, flow_vectors_avgs)
clearmap_io.write(str(dest_path.with_suffix('.tif')), flow_vectors_avgs)
# ############################################################################################################
[docs]
def get_orientation_from_normal_to_surface_local(graph, region_ids, normal_vector_method='mean',
coordinates_name='coordinates'):
"""
Orientation in the cortex relative to the surface normal. (i.e. we fit a plane to the top vertices
and the bottom vertices and take the normal to this plane as the top-to-bottom direction). Alternatively,
we can use the average method which takes the mean of the top and bottom vertices.
Parameters
----------
graph : GraphGt.Graph
The graph to get the orientation from
region_ids : list
The list of the ids of the regions to consider.
normal_vector_method : str
Determines the way to compute the normal vector to the plane best fitting the vertices.
Can be either 'svd' or 'mean'.
If set to svd, we use Singular Value Decomposition (SVD) to find the normal vector to the plane best fitting
the top vertices.
If set to mea, we take the average of the top and bottom vertices.
coordinates_name : str
The coordinates space to use. Typically, 'coordinates', 'coordinates_atlas',
'coordinates_scaled' or 'coordinates_MRI'
Returns
-------
radial_orientations, planar_orientations, reference_norms, edge_lengths
(all ndarrays)
"""
radial_orientations = np.zeros(graph.n_edges)
planar_orientations = np.zeros(graph.n_edges)
label = graph.vertex_annotation()
for id_ in region_ids: # WARNING: could use second element of tuple which is the level
vertex_filter = label == id_
sub_graph = graph.sub_graph(vertex_filter=vertex_filter) # PERFORMANCE: check if we could use view instead of sub_graph
r, p = get_orientation_from_normal_to_surface_global(sub_graph, sub_graph, # FIXME: check if sub_graph is correct
normal_vector_method=normal_vector_method,
coordinates_name=coordinates_name)
edge_filter = vertex_to_edge_property(graph, vertex_filter)
radial_orientations[edge_filter] = r
planar_orientations[edge_filter] = p
return abs(radial_orientations), abs(planar_orientations)
[docs]
def get_orientation_from_normal_to_surface_global(graph, ref_graph, normal_vector_method='mean',
coordinates_name='coordinates'):
"""
Get the orientation of the graph relative to the surface normal.
Parameters
----------
graph : GraphGt.Graph
The graph to get the orientation from
ref_graph
normal_vector_method : str
The method to use to compute the normal vector to the plane best fitting the vertices.
Can be either 'svd' or 'mean'.
coordinates_name : str
The coordinates space to use. Typically, 'coordinates', 'coordinates_atlas',
'coordinates_scaled' or 'coordinates_MRI'
Returns
-------
"""
vertex_coordinates = get_vertex_coordinates(graph, coordinates_name=coordinates_name)
normed_edge_vect = coordinates_to_edge_vectors(vertex_coordinates, graph.edge_connectivity(), normalized=True)
top_to_bottom_dist = get_top_to_bottom_dist(ref_graph, normal_vector_method=normal_vector_method)
radial_orientations = np.dot(top_to_bottom_dist, normed_edge_vect)
planar_orientations = np.sqrt(1 - radial_orientations ** 2)
# edge_vect = coordinates_to_edge_vectors(vertex_coordinates, connectivity, normalized=False)
# reference_norms = np.linalg.norm(edge_vect[edge_vect.shape[0] - 1])
return abs(radial_orientations), abs(planar_orientations)
[docs]
def get_top_to_bottom_dist(graph, normal_vector_method='svd', coordinates_name='coordinates'):
"""
Get the top to bottom distance of the graph.
This can be computed using the minimum distance method or by taking the mean of the top and bottom vertices.
Parameters
----------
graph : GraphGt.Graph
The graph to get the top to bottom distance from
normal_vector_method : str
The method to use to compute the top to bottom distance.
Can be either 'svd' or 'mean'.
If set to 'svd', we identify the vertices in the graph that are closer to the surface than a certain threshold
(min_dist, which is the minimum distance to the surface plus 1.5).
We then calculate the centroid of these vertices and uses Singular Value Decomposition (SVD)
to find the normal vector to the plane best fitting these vertices.
This normal vector is considered as the top-to-bottom direction.
If set to 'mean', the function computes the top-to-bottom distance by taking the mean of the coordinates
of the top and bottom vertices in the graph.
The top vertices are those with a distance to the surface less than or equal to the minimum distance plus 1,
and the bottom vertices are those with a distance to the surface greater than or equal
to the maximum distance minus 1.
The top-to-bottom direction is then the normalized vector from the mean of the top vertices
to the mean of the bottom vertices.
coordinates_name : str
The coordinates space to use. Typically, 'coordinates', 'coordinates_atlas',
'coordinates_scaled' or 'coordinates_MRI'
Returns
-------
"""
dist = graph.vertex_property('distance_to_surface')
coordinates = get_vertex_coordinates(graph, coordinates_name)
if normal_vector_method == 'svd':
margin = 1.5
top_vertices = dist < np.min(dist) + margin
coords = coordinates[top_vertices]
centroid = coords.mean(axis=0)
# run SVD
_, _, vh = np.linalg.svd(coords - centroid)
# unitary normal vector
top_to_bottom_dist = vh[2, :] # FIXME: why 2 ?
elif normal_vector_method == 'mean':
margin = 1 # FIXME: different margin
top_vertices = dist <= np.min(dist) + margin
bottom_vertices = dist >= np.max(dist) - margin
top_to_bottom_dist = (np.mean(coordinates[bottom_vertices], axis=0) -
np.mean(coordinates[top_vertices], axis=0))
top_to_bottom_dist /= np.linalg.norm(top_to_bottom_dist) # preprocessing.normalize(top_to_bottom_dist, norm='l2')
else:
raise ValueError(f'Unknown method {normal_vector_method}, expected "svd" or "mean"')
return top_to_bottom_dist.T
# FIXME: extract part that modifies the graph to add artery and vein
[docs]
def generalized_radial_planar_orientation(graph, work_dir, min_radius, reference_samples, corrected=True,
distance_to_surface=True,
mode='arteryvein', average=False, dvlpmt=True, min_dist_to_surface=7,
coordinates_name='coordinates', orientation_criterion='distance'):
"""
Use Flow interpolation to compute the orientation of the vessels in the graph.
This is based on the local orientation of the arteries
Parameters
----------
graph
work_dir
min_radius
reference_samples
corrected
distance_to_surface
mode
average
dvlpmt
min_dist_to_surface
coordinates_name
orientation_criterion
Returns
-------
"""
artery = edge_to_vertex_property(graph, 'artery')
vein = edge_to_vertex_property(graph, 'vein')
distance_to_surface = graph.vertex_property('distance_to_surface') if distance_to_surface else None
if average:
reference_arteries_edge_coords, reference_arteries_edge_vectors = get_reference_streamlines(work_dir, reference_samples, mode) # FIXME: this could be computed in parent function and passed as argument
capillaries_edge_coordinates, capillaries_edge_vectors, filter = (
get_artery_vein_edge_coords_and_vectors(graph, mode, artery, vein,
min_radius=min_radius,
distance_to_surface=distance_to_surface, min_dist_to_surface=min_dist_to_surface,
coordinates_name=coordinates_name, orientation_criterion=orientation_criterion,
invert=True, return_filter=True))
tmp_angles = [] # The angles to all the reference samples
# PERFORMANCE: should use multiple cores
# For each reference sample, interpolate its arteries vector onto the edges of the tested graph and compute the angle
for arteries_edge_coordinates, arteries_graph_vector in zip(reference_arteries_edge_coords, reference_arteries_edge_vectors):
flow_vectors = interpolate_vectors(arteries_graph_vector, arteries_edge_coordinates, capillaries_edge_coordinates)
tmp_angles.append(angle_between_vectors(flow_vectors, capillaries_edge_vectors))
sub_angles = np.nanmedian(np.array(tmp_angles), axis=0)
angles = np.full(graph.n_edges, np.nan)
indices = np.where(filter)[0]
angles[indices] = sub_angles # Cast back to the original graph shape
else:
edge_coordinates, edge_vectors = graph_to_edge_coords_and_vectors(graph, coordinates_name=coordinates_name, # FIXME: is this capillaries ? do we need to invert ?
orientation_criterion=orientation_criterion)
arteries_edge_coordinates, arteries_edge_vectors = (
get_artery_vein_edge_coords_and_vectors(graph, mode, artery, vein,
min_radius=min_radius, clean_graph=(corrected and not dvlpmt),
distance_to_surface=distance_to_surface, min_dist_to_surface=min_dist_to_surface,
coordinates_name=coordinates_name, orientation_criterion=orientation_criterion))
flow_vectors = interpolate_vectors(arteries_edge_vectors, arteries_edge_coordinates, edge_coordinates)
angles = angle_between_vectors(flow_vectors, edge_vectors)
return angles, graph # FIXME: why return graph ? It shouldn't be modified
