# -*- coding: utf-8 -*-
# -*- python -*-
#
# PropertyTopomesh
#
# Copyright 2014-2016 INRIA - CIRAD - INRA
#
# File author(s): Guillaume Cerutti <guillaume.cerutti@inria.fr>
#
# File contributor(s): Guillaume Cerutti <guillaume.cerutti@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/
#
###############################################################################
from time import time
import numpy as np
from scipy.cluster.vq import vq
try:
import vtk
except ImportError:
print("VTK needs to be installed to use these functionalities")
raise
from copy import deepcopy
from cellcomplex.utils import array_dict
from cellcomplex.property_topomesh import PropertyTopomesh
from cellcomplex.triangular_mesh import TriangularMesh
from cellcomplex.property_topomesh.analysis import compute_topomesh_property, is_triangular
from cellcomplex.property_topomesh.extraction import star_interface_topomesh
triangle_edge_list = np.array([[1, 2],[0, 2],[0, 1]])
[docs]def topomesh_to_triangular_mesh(input_topomesh, degree=3, wids=None, coef=1.0, mesh_center=None, epidermis=False, cell_edges=False, property_name=None, property_degree=None):
topomesh = deepcopy(input_topomesh)
if not is_triangular(topomesh) and (degree > 1):
topomesh = star_interface_topomesh(topomesh)
start_time = time()
print("--> Creating triangular mesh")
triangular_mesh = TriangularMesh()
if property_name is not None:
if property_degree is None:
property_degree = degree
try:
if not topomesh.has_wisp_property(property_name, degree=property_degree, is_computed=True):
compute_topomesh_property(topomesh, property_name, degree=property_degree)
assert len(topomesh.wisp_property(property_name, property_degree).keys()) == topomesh.nb_wisps(property_degree)
except:
property_name = None
if mesh_center is None:
mesh_center = topomesh.wisp_property('barycenter', 0).values().mean(axis=0)
else:
mesh_center = np.array(mesh_center)
considered_start_time = time()
if wids is None:
wids = list(topomesh.wisps(degree))
considered_elements = {}
for d in range(4):
if d == degree:
considered_elements[d] = list(wids)
elif d < degree:
considered_elements[d] = list(np.unique(np.concatenate([list(topomesh.borders(degree, w, degree - d)) for w in wids])))
else:
considered_elements[d] = list(np.unique(np.concatenate([list(topomesh.regions(degree, w, d - degree)) for w in wids])))
if degree > 1:
compute_topomesh_property(topomesh, 'vertices', degree=2)
compute_topomesh_property(topomesh, 'faces', degree=3)
compute_topomesh_property(topomesh, 'cells', degree=2)
cell_triangles = np.concatenate(topomesh.wisp_property('faces', 3).values(considered_elements[3])).astype(int)
# cell_triangles = [t for t in cell_triangles if t in considered_elements[2]]
considered_end_time = time()
print(" --> Filtering out mesh elements [", considered_end_time - considered_start_time, "s]")
if degree == 3:
if property_name is not None:
property_data = topomesh.wisp_property(property_name, property_degree).values(considered_elements[3])
else:
property_data = np.array(considered_elements[3])
vertices_positions = []
triangle_vertices = []
triangle_topomesh_cells = []
vertices_topomesh_vertices = []
# triangle_topomesh_triangles = []
if property_data.ndim == 1 or property_degree < 3:
for c in considered_elements[3]:
if len(list(topomesh.borders(3, c, 3))) > 0:
cell_center = topomesh.wisp_property('barycenter', 0).values(list(topomesh.borders(3, c, 3))).mean(axis=0)
cell_vertices_position = cell_center + coef * (topomesh.wisp_property('barycenter', 0).values(list(topomesh.borders(3, c, 3))) - cell_center) - mesh_center
cell_vertices_index = array_dict(len(vertices_positions) + np.arange(len(list(topomesh.borders(3, c, 3)))), list(topomesh.borders(3, c, 3)))
vertices_positions += list(cell_vertices_position)
vertices_topomesh_vertices += list(topomesh.borders(3, c, 3))
triangle_vertices += list(cell_vertices_index.values(topomesh.wisp_property('vertices', 2).values(topomesh.wisp_property('faces', 3)[c])))
triangle_topomesh_cells += list(c * np.ones_like(topomesh.wisp_property('faces', 3)[c]))
# triangle_topomesh_triangles += topomesh.wisp_property('faces',3)[c]
print(len(cell_triangles), "Cell triangles")
print(len(triangle_vertices), "Triangle vertices")
vertices_positions = array_dict(vertices_positions, np.arange(len(vertices_positions)))
vertices_topomesh_vertices = array_dict(vertices_topomesh_vertices, np.arange(len(vertices_positions)))
if epidermis:
compute_topomesh_property(topomesh, 'epidermis', 2)
epidermis_triangles = topomesh.wisp_property('epidermis', 2).values(cell_triangles)
triangle_vertices = array_dict(np.array(triangle_vertices)[epidermis_triangles], np.arange(len(cell_triangles[epidermis_triangles])))
triangle_topomesh_cells = array_dict(np.array(triangle_topomesh_cells)[epidermis_triangles], np.arange(len(cell_triangles[epidermis_triangles])))
triangle_topomesh_triangles = array_dict(cell_triangles[epidermis_triangles], np.arange(len(cell_triangles[epidermis_triangles])))
else:
triangle_vertices = array_dict(triangle_vertices, np.arange(len(cell_triangles)))
triangle_topomesh_cells = array_dict(triangle_topomesh_cells, np.arange(len(cell_triangles)))
triangle_topomesh_triangles = array_dict(cell_triangles, np.arange(len(cell_triangles)))
edge_topomesh_edges = {}
triangular_mesh.points = vertices_positions.to_dict()
triangular_mesh.triangles = triangle_vertices.to_dict()
if property_name is not None:
if property_degree == 2:
triangle_topomesh_triangle_property = array_dict(topomesh.wisp_property(property_name, property_degree).values(triangle_topomesh_triangles.values()),
triangle_topomesh_triangles.keys())
triangular_mesh.triangle_data = triangle_topomesh_triangle_property.to_dict()
elif property_degree == 0:
vertex_topomesh_vertex_property = array_dict(topomesh.wisp_property(property_name, property_degree).values(vertices_topomesh_vertices.values()), vertices_topomesh_vertices.keys())
triangular_mesh.point_data = vertex_topomesh_vertex_property.to_dict()
elif property_degree == 3:
triangle_topomesh_cell_property = array_dict(topomesh.wisp_property(property_name, property_degree).values(triangle_topomesh_cells.values()), triangle_topomesh_cells.keys())
triangular_mesh.triangle_data = triangle_topomesh_cell_property.to_dict()
else:
triangular_mesh.triangle_data = triangle_topomesh_cells.to_dict()
else:
for c in considered_elements[3]:
if len(list(topomesh.borders(3, c, 3))) > 0:
cell_center = topomesh.wisp_property('barycenter', 0).values(list(topomesh.borders(3, c, 3))).mean(axis=0)
vertices_positions += [cell_center]
vertices_positions = array_dict(vertices_positions, np.array([c for c in topomesh.wisps(3) if len(list(topomesh.borders(3, c, 3))) > 0]))
vertices_topomesh_vertices = {}
edge_topomesh_edges = {}
triangle_topomesh_triangles = {}
triangle_topomesh_cells = {}
cell_property = array_dict(topomesh.wisp_property(property_name, property_degree).values(vertices_positions.keys()), vertices_positions.keys())
triangular_mesh.points = vertices_positions.to_dict()
triangular_mesh.point_data = cell_property
triangular_mesh.triangles = {}
elif degree == 2:
vertices_positions = []
triangle_vertices = []
vertices_topomesh_vertices = []
for t in cell_triangles:
triangle_center = topomesh.wisp_property('barycenter', 0).values(list(topomesh.borders(2, t, 2))).mean(axis=0)
triangle_vertices_position = triangle_center + coef * (topomesh.wisp_property('barycenter', 0).values(list(topomesh.borders(2, t, 2))) - triangle_center) - mesh_center
triangle_vertices_index = array_dict(len(vertices_positions) + np.arange(3), list(topomesh.borders(2, t, 2)))
vertices_positions += list(triangle_vertices_position)
vertices_topomesh_vertices += list(topomesh.borders(2, t, 2))
triangle_vertices += list([triangle_vertices_index.values(topomesh.wisp_property('vertices', 2)[t])])
vertices_positions = array_dict(vertices_positions, np.arange(len(vertices_positions)))
vertices_topomesh_vertices = array_dict(vertices_topomesh_vertices, np.arange(len(vertices_positions)))
triangle_topomesh_cells = np.concatenate([c * np.ones_like([t for t in topomesh.wisp_property('faces', 3)[c] if t in considered_elements[2]]) for c in considered_elements[3]]).astype(int)
if epidermis:
compute_topomesh_property(topomesh, 'epidermis', 2)
epidermis_triangles = topomesh.wisp_property('epidermis', 2).values(cell_triangles)
triangle_vertices = array_dict(np.array(triangle_vertices)[epidermis_triangles], np.arange(len(cell_triangles[epidermis_triangles])))
triangle_topomesh_cells = array_dict(np.array(triangle_topomesh_cells)[epidermis_triangles], np.arange(len(cell_triangles[epidermis_triangles])))
triangle_topomesh_triangles = array_dict(cell_triangles[epidermis_triangles], np.arange(len(cell_triangles[epidermis_triangles])))
else:
triangle_vertices = array_dict(triangle_vertices, np.arange(len(cell_triangles)))
triangle_topomesh_cells = array_dict(triangle_topomesh_cells, np.arange(len(cell_triangles)))
triangle_topomesh_triangles = array_dict(cell_triangles, np.arange(len(cell_triangles)))
edge_topomesh_edges = {}
triangular_mesh.points = vertices_positions.to_dict()
triangular_mesh.triangles = triangle_vertices.to_dict()
if property_name is not None:
if property_degree == 2:
triangle_topomesh_triangle_property = array_dict(topomesh.wisp_property(property_name, property_degree).values(triangle_topomesh_triangles.values()),
triangle_topomesh_triangles.keys())
triangular_mesh.triangle_data = triangle_topomesh_triangle_property.to_dict()
elif property_degree == 0:
vertex_topomesh_vertex_property = array_dict(topomesh.wisp_property(property_name, property_degree).values(vertices_topomesh_vertices.values()), vertices_topomesh_vertices.keys())
triangular_mesh.point_data = vertex_topomesh_vertex_property.to_dict()
elif property_degree == 3:
triangle_topomesh_cell_property = array_dict(topomesh.wisp_property(property_name, property_degree).values(triangle_topomesh_cells.values()), triangle_topomesh_cells.keys())
triangular_mesh.triangle_data = triangle_topomesh_cell_property.to_dict()
else:
triangular_mesh.triangle_data = triangle_topomesh_cells.to_dict()
elif degree == 1:
compute_topomesh_property(topomesh, 'vertices', degree=1)
vertices_positions = array_dict(topomesh.wisp_property('barycenter', 0).values(considered_elements[0]), considered_elements[0])
edge_vertices = array_dict(topomesh.wisp_property('vertices', 1).values(considered_elements[1]), considered_elements[1])
triangular_mesh.points = vertices_positions.to_dict()
triangular_mesh.edges = edge_vertices.to_dict()
if property_name is not None:
if property_degree == 1:
edge_property = array_dict(topomesh.wisp_property(property_name, property_degree).values(considered_elements[1]), considered_elements[1])
triangular_mesh.edge_data = edge_property.to_dict()
triangle_topomesh_cells = {}
triangle_topomesh_triangles = {}
edge_topomesh_edges = dict(list(zip(triangular_mesh.edges.keys(), triangular_mesh.edges.keys())))
vertices_topomesh_vertices = {}
if cell_edges:
compute_topomesh_property(topomesh, 'epidermis', 1)
compute_topomesh_property(topomesh, 'cells', 1)
edge_n_cells = np.array(list(map(len, topomesh.wisp_property('cells', 1).values())))
edge_is_cell_edge = edge_n_cells > 2
edge_is_cell_edge = edge_is_cell_edge | (edge_n_cells > 1) & (topomesh.wisp_property('epidermis', 1).values())
edge_is_cell_edge = array_dict(edge_is_cell_edge, list(topomesh.wisps(1)))
edges_to_remove = []
for eid in triangular_mesh.edges:
if not edge_is_cell_edge[eid]:
edges_to_remove += [eid]
for eid in edges_to_remove:
del triangular_mesh.edges[eid]
del edge_topomesh_edges[eid]
if eid in triangular_mesh.edge_data:
del triangular_mesh.edge_data[eid]
elif degree == 0:
vertices_positions = array_dict(topomesh.wisp_property('barycenter', 0).values(considered_elements[0]), considered_elements[0])
triangular_mesh.points = vertices_positions.to_dict()
if property_name is not None:
if property_degree == 0:
vertex_property = array_dict(topomesh.wisp_property(property_name, property_degree).values(considered_elements[0]), considered_elements[0])
triangular_mesh.point_data = vertex_property.to_dict()
triangle_topomesh_cells = {}
triangle_topomesh_triangles = {}
edge_topomesh_edges = {}
vertices_topomesh_vertices = dict(list(zip(triangular_mesh.points.keys(), triangular_mesh.points.keys())))
mesh_element_matching = {}
mesh_element_matching[0] = vertices_topomesh_vertices
mesh_element_matching[1] = edge_topomesh_edges
mesh_element_matching[2] = triangle_topomesh_triangles
mesh_element_matching[3] = triangle_topomesh_cells
end_time = time()
print("<-- Creating triangular mesh [", end_time - start_time, "s]")
return triangular_mesh, mesh_element_matching
[docs]def vtk_polydata_to_triangular_mesh(polydata):
mesh = TriangularMesh()
start_time = time()
print("--> Creating vertices")
polydata_point_data = None
if polydata.GetPointData().GetNumberOfComponents() > 0:
polydata_point_data = polydata.GetPointData().GetArray(0)
for v in range(polydata.GetPoints().GetNumberOfPoints()):
mesh.points[v] = polydata.GetPoints().GetPoint(v)
if polydata_point_data is not None:
mesh.point_data[v] = polydata_point_data.GetTuple(v)
end_time = time()
print("<-- Creating vertices [",end_time-start_time,"s]")
start_time = time()
print("--> Creating triangles")
triangles = {}
triangle_data = {}
polydata_cell_data = None
if polydata.GetCellData().GetNumberOfComponents() > 0:
polydata_cell_data = polydata.GetCellData().GetArray(0)
for t in range(polydata.GetNumberOfCells()):
mesh.triangles[t] = np.sort([polydata.GetCell(t).GetPointIds().GetId(i) for i in range(3)])
if polydata_cell_data is not None:
mesh.triangle_data[t] = polydata_cell_data.GetTuple(t)
end_time = time()
print("<-- Creating triangles [",end_time-start_time,"s]")
return mesh
[docs]def vtk_polydata_to_cell_triangular_meshes(polydata):
mesh = {}
polydata_cell_data = polydata.GetCellData().GetArray(0)
triangle_cell_start_time = time()
print(" --> Listing triangles")
print(" - ", polydata.GetNumberOfCells(), " triangles")
# polydata_triangles = np.sort([[polydata.GetCell(t).GetPointIds().GetId(i) for i in range(3)] for t in range(polydata.GetNumberOfCells())])
polydata_triangles = np.array([[polydata.GetCell(t).GetPointIds().GetId(i) for i in range(3)] for t in range(polydata.GetNumberOfCells())])
triangle_cell_end_time = time()
print(" <-- Listing triangles [", triangle_cell_end_time - triangle_cell_start_time, "s]")
triangle_cell_start_time = time()
print(" --> Listing triangle cells")
triangle_cell = np.array([polydata_cell_data.GetTuple(t)[0] for t in range(polydata.GetNumberOfCells())], np.uint16)
triangle_cell_end_time = time()
print(" <-- Listing triangle cells [", triangle_cell_end_time - triangle_cell_start_time, "s]")
start_time = time()
print(" --> Creating cell meshes")
for c in np.unique(triangle_cell):
mesh[c] = TriangularMesh()
cell_triangles = np.arange(polydata.GetNumberOfCells())[np.where(triangle_cell == c)]
# cell_triangle_points = np.sort([[polydata.GetCell(t).GetPointIds().GetId(i) for i in range(3)] for t in cell_triangles])
cell_triangle_points = np.array([[polydata.GetCell(t).GetPointIds().GetId(i) for i in range(3)] for t in cell_triangles])
cell_vertices = np.sort(np.unique(cell_triangle_points))
mesh[c].points = array_dict(np.array([polydata.GetPoints().GetPoint(v) for v in cell_vertices]), cell_vertices).to_dict()
mesh[c].triangles = array_dict(cell_triangle_points, cell_triangles).to_dict()
mesh[c].triangle_data = array_dict(np.ones_like(cell_triangles) * c, cell_triangles).to_dict()
end_time = time()
print(" <-- Creating cell meshes [", end_time - start_time, "s]")
return mesh
[docs]def vtk_polydata_to_topomesh(polydata):
topomesh = PropertyTopomesh(3)
start_time = time()
print("--> Creating vertices")
# polydata_points = np.array([polydata.GetPoints().GetPoint(i) for i in range(polydata.GetPoints().GetNumberOfPoints())])
# unique_points = array_unique(polydata_points)
# n_points = unique_points.shape[0]
# polydata_point_point_matching = (vq(polydata_points,unique_points)[0])
vertex_positions = {}
for v in range(polydata.GetPoints().GetNumberOfPoints()):
# for v in range(n_points):
vertex_positions[v] = polydata.GetPoints().GetPoint(v)
# vertex_positions[v] = unique_points[v]
topomesh.add_wisp(0, v)
end_time = time()
print("<-- Creating vertices [", end_time - start_time, "s]")
start_time = time()
print("--> Creating edges")
polydata_triangles = np.sort([[polydata.GetCell(t).GetPointIds().GetId(i) for i in range(3)] for t in range(polydata.GetNumberOfCells())])
# polydata_triangles = np.sort(polydata_point_point_matching[polydata_triangles])
polydata_edges = np.unique(np.concatenate(polydata_triangles[:, triangle_edge_list]), axis=0)
for e, edge_vertices in enumerate(polydata_edges):
topomesh.add_wisp(1, e)
for v in edge_vertices:
topomesh.link(1, e, v)
end_time = time()
print("<-- Creating edges [", end_time - start_time, "s]")
start_time = time()
print("--> Linking triangle to edges")
unique_triangles = np.unique(polydata_triangles, axis=0)
n_triangles = unique_triangles.shape[0]
polydata_triangle_triangle_matching = (vq(polydata_triangles, unique_triangles)[0])
polydata_triangle_edge_matching = (vq(np.concatenate(unique_triangles[:, triangle_edge_list]), polydata_edges)[0]).reshape((n_triangles, 3))
end_time = time()
print("<-- Linking triangle to edges[", end_time - start_time, "s]")
start_time = time()
print("--> Creating triangles")
for t in range(n_triangles):
topomesh.add_wisp(2, t)
for e in polydata_triangle_edge_matching[t]:
topomesh.link(2, t, e)
end_time = time()
print("<-- Creating triangles [", end_time - start_time, "s]")
start_time = time()
print("--> Creating cells")
polydata_cell_data = polydata.GetCellData().GetArray(0)
polydata_cells = np.array(np.unique(np.concatenate([polydata_cell_data.GetTuple(t) for t in range(polydata.GetNumberOfCells())])), np.uint16)
for c in polydata_cells:
if c != 1:
topomesh.add_wisp(3, c)
for t in range(polydata.GetNumberOfCells()):
for c in polydata_cell_data.GetTuple(t):
if c != 1:
topomesh.link(3, int(c), polydata_triangle_triangle_matching[t])
end_time = time()
print("<-- Creating cells [", end_time - start_time, "s]")
topomesh.update_wisp_property('barycenter', 0, vertex_positions)
return topomesh
[docs]def topomesh_to_vtk_polydata(topomesh, degree=2, positions=None, topomesh_center=None, coef=1):
if positions is None:
positions = topomesh.wisp_property('barycenter', 0)
if topomesh_center is None:
topomesh_center = np.mean(positions.values(), axis=0)
# topomesh_center = np.array([0,0,0])
print(topomesh_center)
vtk_mesh = vtk.vtkPolyData()
vtk_points = vtk.vtkPoints()
vtk_edges = vtk.vtkCellArray()
vtk_triangles = vtk.vtkCellArray()
vtk_cells = vtk.vtkLongArray()
start_time = time()
print("--> Creating VTK PolyData")
if degree == 3:
for c in topomesh.wisps(3):
cell_points = []
cell_center = np.mean([positions[v] for v in topomesh.borders(3, c, 3)], axis=0)
for v in topomesh.borders(3, c, 3):
position = cell_center + coef * (positions[v] - cell_center) - topomesh_center
position[np.where(np.abs(position) < .001)] = 0.
p = vtk_points.InsertNextPoint(position)
cell_points.append(p)
cell_points = array_dict(cell_points, list(topomesh.borders(3, c, 3)))
for t in topomesh.borders(3, c):
poly = vtk_triangles.InsertNextCell(3)
for v in topomesh.borders(2, t, 2):
vtk_triangles.InsertCellPoint(cell_points[v])
vtk_cells.InsertValue(poly, c)
elif degree == 2:
for t in topomesh.wisps(2):
triangle_points = []
triangle_center = np.mean([positions[v] for v in topomesh.borders(2, t, 2)], axis=0)
for v in topomesh.borders(2, t, 2):
position = triangle_center + coef * (positions[v] - triangle_center) - topomesh_center
position[np.where(np.abs(position) < .001)] = 0.
p = vtk_points.InsertNextPoint(position)
triangle_points.append(p)
triangle_points = array_dict(triangle_points, list(topomesh.borders(2, t, 2)))
poly = vtk_triangles.InsertNextCell(3)
for v in topomesh.borders(2, t, 2):
vtk_triangles.InsertCellPoint(triangle_points[v])
vtk_cells.InsertValue(poly, list(topomesh.regions(2, t))[0])
elif degree == 1:
points = []
for v in topomesh.wisps(0):
position = positions[v]
position[np.where(np.abs(position) < .001)] = 0.
p = vtk_points.InsertNextPoint(position)
points.append(p)
points = array_dict(points, list(topomesh.wisps(0)))
for e in topomesh.wisps(1):
# if topomesh.wisp_property('epidermis',1)[e]:
# if True:
if len(list(topomesh.regions(1, e))) > 2:
c = vtk_edges.InsertNextCell(2)
for v in topomesh.borders(1, e):
vtk_edges.InsertCellPoint(points[v])
print(" --> Linking Mesh")
vtk_mesh.SetPoints(vtk_points)
vtk_mesh.SetPolys(vtk_triangles)
vtk_mesh.SetLines(vtk_edges)
vtk_mesh.GetCellData().SetScalars(vtk_cells)
end_time = time()
print("<-- Creating VTK PolyData [", end_time - start_time, "s]")
return vtk_mesh
[docs]def composed_triangular_mesh(triangular_mesh_dict):
start_time = time()
print("--> Composing triangular mesh...")
mesh = TriangularMesh()
triangle_cell_matching = {}
mesh_points = np.concatenate([triangular_mesh_dict[c].points.keys() for c in triangular_mesh_dict.keys()])
mesh_point_positions = np.concatenate([triangular_mesh_dict[c].points.values() for c in triangular_mesh_dict.keys()])
mesh.points = dict(list(zip(mesh_points,mesh_point_positions)))
mesh_triangles = np.concatenate([triangular_mesh_dict[c].triangles.values() for c in triangular_mesh_dict.keys()])
mesh.triangles = dict(list(zip(np.arange(len(mesh_triangles)),mesh_triangles)))
mesh_cells = np.concatenate([c*np.ones_like(triangular_mesh_dict[c].triangles.keys()) for c in triangular_mesh_dict.keys()])
triangle_cell_matching = dict(list(zip(np.arange(len(mesh_triangles)),mesh_cells)))
end_time = time()
print("<-- Composing triangular mesh [",end_time-start_time,"]")
return mesh, triangle_cell_matching
