Source code for cellcomplex.property_topomesh.utils.tissue_analysis_tools
# -*- coding: utf-8 -*-
# -*- python -*-
#
# PropertyTopomesh
#
# Copyright 2015-2016 INRIA - CIRAD - INRA
#
# File author(s): Guillaume Cerutti <guillaume.cerutti@inria.fr>
#
# File contributor(s): Guillaume Cerutti <guillaume.cerutti@inria.fr>,
# Jonathan Legrand <jonathan.legrand@inria.fr>
#
# Distributed under the Cecill-C License.
# See accompanying file LICENSE.txt or copy at
# http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
#
# OpenaleaLab Website : http://virtualplants.github.io/
#
###############################################################################
import logging
import numpy as np
from time import time as current_time
from cellcomplex.property_topomesh.analysis import compute_topomesh_property
from scipy.cluster.vq import vq
from time import time
[docs]def cell_topological_elements_extraction(img, dimensions=[2,1,0]):
start_time = time()
logging.info("--> Extracting cell topological elements")
shape = np.array(img.shape)
neighborhood_img = []
for x in np.arange(-1,1):
for y in np.arange(-1,1):
for z in np.arange(-1,1):
neighborhood_img.append(img[1+x:shape[0]+x,1+y:shape[1]+y,1+z:shape[2]+z])
neighborhood_img = np.sort(np.transpose(neighborhood_img,(1,2,3,0))).reshape((shape-1).prod(),8)
vertex_coords = np.transpose(np.mgrid[0:shape[0]-1,0:shape[1]-1,0+z:shape[2]-1],(1,2,3,0)).reshape((shape-1).prod(),3)+0.5
non_flat = np.sum(neighborhood_img == np.tile(neighborhood_img[:,:1],(1,8)),axis=1)!=8
neighborhood_img = neighborhood_img[non_flat]
vertex_coords = vertex_coords[non_flat]
neighborhoods = np.array(list(map(np.unique,neighborhood_img)))
neighborhood_size = np.array(list(map(len,neighborhoods)))
neighborhoods = np.array(neighborhoods)
topological_elements = {}
for dimension in dimensions:
element_points = vertex_coords[neighborhood_size==4-dimension]
element_cells = np.array([p for p in neighborhoods[neighborhood_size==4-dimension]],int)
if (dimension==0) & ((neighborhood_size>=5).sum() > 0):
clique_vertex_points = np.concatenate([(p,p) for p in vertex_coords[neighborhood_size==5]])
clique_vertex_cells = np.concatenate([[p[:4],np.concatenate([[p[0]],p[2:]])] for p in neighborhoods[neighborhood_size==5]]).astype(int)
element_points = np.concatenate([element_points,clique_vertex_points])
element_cells = np.concatenate([element_cells,clique_vertex_cells])
unique_cell_elements = np.unique(element_cells,axis=0)
element_matching = vq(element_cells,unique_cell_elements)[0]
topological_elements[dimension] = dict(list(zip([tuple(e) for e in unique_cell_elements],[element_points[element_matching == e] for e,_ in enumerate(unique_cell_elements)])))
end_time = time()
logging.info("<-- Extracting cell topological elements [",end_time - start_time,"s]")
return topological_elements
[docs]def cell_vertex_extraction(img,physical_units=False):
start_time = time()
print("--> Extracting cell vertices")
shape = np.array(img.shape)
if img.ndim==3 and shape[2]>1:
neighborhood_img = []
for x in np.arange(-1,1):
for y in np.arange(-1,1):
for z in np.arange(-1,1):
neighborhood_img.append(img[1+x:shape[0]+x,1+y:shape[1]+y,1+z:shape[2]+z])
neighborhood_img = np.sort(np.transpose(neighborhood_img,(1,2,3,0))).reshape((shape-1).prod(),8)
vertex_coords = np.transpose(np.mgrid[0:shape[0]-1,0:shape[1]-1,0+z:shape[2]-1],(1,2,3,0)).reshape((shape-1).prod(),3)+0.5
non_flat = np.sum(neighborhood_img == np.tile(neighborhood_img[:,:1],(1,8)),axis=1)!=8
neighborhood_img = neighborhood_img[non_flat]
vertex_coords = vertex_coords[non_flat]
neighborhoods = np.array(list(map(np.unique,neighborhood_img)))
neighborhood_size = np.array(list(map(len,neighborhoods)))
neighborhoods = np.array(neighborhoods)
# vertex_coords = np.where(neighborhood_size==4)
#vertex_coords = np.where(neighborhood_size >= 4)
# vertex_points = np.transpose(vertex_coords)+0.5
vertex_points = vertex_coords[neighborhood_size==4]
# vertex_cells = np.array([p for p in neighborhoods[vertex_coords]],int)
vertex_cells = np.array([p for p in neighborhoods[neighborhood_size==4]],int)
#vertex_cells = np.array([p[:4] for p in neighborhoods[vertex_coords]],int)
if (neighborhood_size>=5).sum() > 0:
clique_vertex_points = np.concatenate([(p,p) for p in vertex_coords[neighborhood_size==5]])
clique_vertex_cells = np.concatenate([[p[:4],np.concatenate([[p[0]],p[2:]])] for p in neighborhoods[neighborhood_size==5]]).astype(int)
# clique_vertex_coords = np.where(neighborhood_size==5)
# clique_vertex_points = np.concatenate([(p,p) for p in np.transpose(clique_vertex_coords)])+0.5
# clique_vertex_cells = np.concatenate([[p[:4],np.concatenate([[p[0]],p[2:]])] for p in neighborhoods[clique_vertex_coords]]).astype(int)
# clique_vertex_points = np.concatenate([(p,p,p,p,p) for p in np.transpose(clique_vertex_coords)])+0.5
# clique_vertex_cells = np.concatenate([[np.concatenate([p[:i],p[i+1:]]) for i in range(5)] for p in neighborhoods[clique_vertex_coords]]).astype(int)
vertex_points = np.concatenate([vertex_points,clique_vertex_points])
vertex_cells = np.concatenate([vertex_cells,clique_vertex_cells])
unique_cell_vertices = np.unique(vertex_cells,axis=0)
vertices_matching = vq(vertex_cells,unique_cell_vertices)[0]
unique_cell_vertex_points = np.array([np.mean(vertex_points[vertices_matching == v],axis=0) for v in range(len(unique_cell_vertices))])
if physical_units:
cell_vertex_dict = dict(zip([tuple(c) for c in unique_cell_vertices],unique_cell_vertex_points*np.array(img.voxelsize)))
else:
cell_vertex_dict = dict(zip([tuple(c) for c in unique_cell_vertices],unique_cell_vertex_points))
elif img.ndim==2 or (img.ndim==3 and shape[2]==1):
if img.ndim==3:
img = img[:,:,0]
shape = np.array(img.shape)
neighborhood_img = []
for x in np.arange(-1,1):
for y in np.arange(-1,1):
neighborhood_img.append(img[1+x:shape[0]+x,1+y:shape[1]+y])
neighborhood_img = np.sort(np.transpose(neighborhood_img,(1,2,0))).reshape((shape-1).prod(),4)
neighborhoods = np.array(list(map(np.unique,neighborhood_img)))
neighborhood_size = np.array(list(map(len,neighborhoods))).reshape(shape[0]-1,shape[1]-1)
neighborhoods = np.array(neighborhoods).reshape(shape[0]-1,shape[1]-1)
vertex_coords = np.where(neighborhood_size==3)
vertex_points = np.transpose(vertex_coords)+0.5
vertex_cells = np.array([p for p in neighborhoods[vertex_coords]],int)
if (neighborhood_size>=4).sum() > 0:
clique_vertex_coords = np.where(neighborhood_size==4)
clique_vertex_points = np.concatenate([(p,p) for p in np.transpose(clique_vertex_coords)])+0.5
clique_vertex_cells = np.concatenate([[p[:3],np.concatenate([[p[0]],p[2:]])] for p in neighborhoods[clique_vertex_coords]]).astype(int)
vertex_points = np.concatenate([vertex_points,clique_vertex_points])
vertex_cells = np.concatenate([vertex_cells,clique_vertex_cells])
unique_cell_vertices = np.unique(vertex_cells,axis=0)
vertices_matching = vq(vertex_cells,unique_cell_vertices)[0]
unique_cell_vertex_points = np.array([np.mean(vertex_points[vertices_matching == v],axis=0) for v in range(len(unique_cell_vertices))])
cell_vertex_dict = dict(list(zip([tuple(c) for c in unique_cell_vertices],list(unique_cell_vertex_points))))
else:
cell_vertex_dict = {}
end_time = time()
print("<-- Extracting cell vertices (",len(cell_vertex_dict),"vertices ) [",end_time - start_time,"s]")
return cell_vertex_dict
[docs]def voxel_cell_vertex_extraction(img,**kwargs):
shape = np.array(img.shape)
neighborhood_img = []
for x in np.arange(-1,2):
for y in np.arange(-1,2):
for z in np.arange(-1,2):
neighborhood_img.append(img[1+x:shape[0]-1+x,1+y:shape[1]-1+y,1+z:shape[2]-1+z])
neighborhood_img = np.sort(np.transpose(neighborhood_img,(1,2,3,0))).reshape((shape-2).prod(),27)
neighborhoods = np.array(list(map(np.unique,neighborhood_img)))
neighborhood_size = np.array(list(map(len,neighborhoods))).reshape(shape[0]-2,shape[1]-2,shape[2]-2)
neighborhoods = np.array(neighborhoods).reshape(shape[0]-2,shape[1]-2,shape[2]-2)
vertex_coords = np.where(neighborhood_size==4)
vertex_points = np.transpose(vertex_coords)+1
vertex_cells = np.array([p for p in neighborhoods[vertex_coords]],int)
unique_cell_vertices = np.unique(vertex_cells,axis=0)
vertices_matching = vq(vertex_cells,unique_cell_vertices)[0]
unique_cell_vertex_points = np.array([np.mean(vertex_points[vertices_matching == v],axis=0) for v in range(len(unique_cell_vertices))])
cell_vertex_dict = dict(list(zip([tuple(c) for c in unique_cell_vertices],list(unique_cell_vertex_points))))
return cell_vertex_dict
[docs]def topomesh_cell_vertex_extraction(topomesh, background_label=1):
compute_topomesh_property(topomesh, 'cells', degree=0)
vertex_cell_neighbours = topomesh.wisp_property('cells', degree=0)
surface_faces = [f for f in topomesh.wisps(2) if topomesh.nb_regions(2, f) == 1]
surface_vertices = [v for v in topomesh.wisps(0) if np.any([f in surface_faces for f in topomesh.regions(0, v, 2)])]
start_time = current_time()
print("--> Computing mesh cell vertices")
mesh_cell_vertex = {}
for v in topomesh.wisps(0):
if len(vertex_cell_neighbours[v]) == 5:
for k in range(5):
vertex_cell_labels = tuple([c for c in vertex_cell_neighbours[v]][:k]) + tuple([c for c in vertex_cell_neighbours[v]][k + 1:])
if vertex_cell_labels not in mesh_cell_vertex:
mesh_cell_vertex[vertex_cell_labels] = v
if len(vertex_cell_neighbours[v]) == 4:
vertex_cell_labels = tuple([c for c in vertex_cell_neighbours[v]])
mesh_cell_vertex[vertex_cell_labels] = v
if v in surface_vertices:
# and count_l1_cells(vertex_cell_neighbours[v]) == 4:
for k in range(4):
vertex_cell_labels = (background_label,) + tuple([c for c in vertex_cell_neighbours[v]][:k]) + tuple([c for c in vertex_cell_neighbours[v]][k + 1:])
if vertex_cell_labels not in mesh_cell_vertex:
mesh_cell_vertex[vertex_cell_labels] = v
if (len(vertex_cell_neighbours[v]) == 3) and (v in surface_vertices):
# and (count_l1_cells(vertex_cell_neighbours[v]) == 3):
vertex_cell_labels = (background_label,) + tuple([c for c in vertex_cell_neighbours[v]])
mesh_cell_vertex[vertex_cell_labels] = v
end_time = current_time()
print("<-- Computing mesh cell vertices [", end_time - start_time, "s]")
return mesh_cell_vertex
[docs]def move_topomesh_cell_vertices(topomesh, cell_vertex_positions, background_label=1, coef=0.5):
mesh_cell_vertex = topomesh_cell_vertex_extraction(topomesh, background_label=background_label)
# cell_vertex_matching = vq(np.sort(array_dict(cell_vertex_positions).keys()),np.sort(array_dict(mesh_cell_vertex).keys()))
cell_vertex_matching = vq(np.sort(np.array(list(cell_vertex_positions.keys()))), np.sort(np.array(list(mesh_cell_vertex.keys()))))
matched_image_index = np.where(cell_vertex_matching[1] == 0)[0]
matched_mesh_index = cell_vertex_matching[0][np.where(cell_vertex_matching[1] == 0)[0]]
matched_cell_vertex_positions = np.array(list(cell_vertex_positions.values()))[matched_image_index]
matched_keys = np.sort(np.array(list(cell_vertex_positions.keys())))[matched_image_index]
matched_mesh_vertices = np.array(list(mesh_cell_vertex.values()))[cell_vertex_matching[0][np.where(cell_vertex_matching[1] == 0)[0]]]
matched_keys = np.sort(np.array(list(mesh_cell_vertex.keys())))[matched_mesh_index]
vertex_positions = topomesh.wisp_property('barycenter', 0).to_dict()
for i, v in enumerate(matched_mesh_vertices):
if not np.isnan(matched_cell_vertex_positions[i]).any():
vertex_positions[v] = coef * matched_cell_vertex_positions[i] + (1 - coef) * vertex_positions[v]
topomesh.update_wisp_property('barycenter', 0, vertex_positions)
