# -*- coding: utf-8 -*-
# -*- python -*-
#
# PropertyTopomesh
#
# Copyright 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/
#
###############################################################################
import numpy as np
from cellcomplex.utils import array_dict
from cellcomplex.property_topomesh.creation import triangle_topomesh, quad_topomesh, poly_topomesh, dual_topomesh
from cellcomplex.property_topomesh.analysis import compute_topomesh_property, compute_topomesh_triangle_properties
from cellcomplex.property_topomesh.extraction import star_interface_topomesh
from cellcomplex.property_topomesh.optimization import topomesh_triangle_split, property_topomesh_vertices_deformation
from cellcomplex.property_topomesh.composition import append_topomesh, fuse_topomesh_identical_vertices
from cellcomplex.property_topomesh.utils.delaunay_tools import delaunay_triangulation
[docs]def square_topomesh(side_length = 1):
points = {}
points[0] = [0,0,0]
points[1] = [side_length,0,0]
points[2] = [0,side_length,0]
points[3] = [side_length,side_length,0]
triangles = [[0,1,2],[3,2,1]]
return triangle_topomesh(triangles, points)
[docs]def hexagon_topomesh(side_length = 1, angle=0.):
points = {}
points[0] = [0,0,0]
for p in range(6):
points[p+1] = [side_length*np.cos(p*np.pi/3.+angle),side_length*np.sin(p*np.pi/3.+angle),0]
triangles = []
for p in range(6):
triangles += [[0,p + 1,(p+1)%6 + 1]]
return triangle_topomesh(triangles, points)
[docs]def icosahedron_topomesh(center=np.array([0,0,0],float),radius=10.):
ico_points = {}
ico_points[0] = center + radius*np.array([ 0. , -1. , 0. ])
ico_points[1] = center + radius*np.array([ 0.7236 , -0.447215, 0.52572 ])
ico_points[2] = center + radius*np.array([-0.276385, -0.447215, 0.85064 ])
ico_points[3] = center + radius*np.array([-0.894425, -0.447215, 0. ])
ico_points[4] = center + radius*np.array([-0.276385, -0.447215, -0.85064 ])
ico_points[5] = center + radius*np.array([ 0.7236 , -0.447215, -0.52572 ])
ico_points[6] = center + radius*np.array([ 0.276385, 0.447215, 0.85064 ])
ico_points[7] = center + radius*np.array([-0.7236 , 0.447215, 0.52572 ])
ico_points[8] = center + radius*np.array([-0.7236 , 0.447215, -0.52572 ])
ico_points[9] = center + radius*np.array([ 0.276385, 0.447215, -0.85064 ])
ico_points[10] = center + radius*np.array([ 0.894425, 0.447215, 0. ])
ico_points[11] = center + radius*np.array([ 0. , 1. , 0. ])
ico_triangles = [[ 0, 1, 2],
[ 0, 1, 5],
[ 0, 2, 3],
[ 0, 3, 4],
[ 0, 4, 5],
[ 1, 5, 10],
[ 1, 2, 6],
[ 2, 3, 7],
[ 3, 4, 8],
[ 4, 5, 9],
[ 1, 6, 10],
[ 2, 6, 7],
[ 3, 7, 8],
[ 4, 8, 9],
[ 5, 9, 10],
[ 6, 10, 11],
[ 6, 7, 11],
[ 7, 8, 11],
[ 8, 9, 11],
[ 9, 10, 11]]
return triangle_topomesh(ico_triangles, ico_points)
[docs]def truncated_octahedron_topomesh(center=np.array([0,0,0],float),radius=10.,ratio=1.):
points = []
points.append(np.array([0,-1,-1]))
points.append(np.array([0,-1,1]))
points.append(np.array([0,1,1]))
points.append(np.array([0,1,-1]))
points.append(np.array([np.sqrt(2),0,0]))
points.append(np.array([-np.sqrt(2),0,0]))
triangles = []
triangles.append(np.array([0,1,4]))
triangles.append(np.array([1,2,4]))
triangles.append(np.array([2,3,4]))
triangles.append(np.array([3,0,4]))
triangles.append(np.array([1,0,5]))
triangles.append(np.array([2,1,5]))
triangles.append(np.array([3,2,5]))
triangles.append(np.array([0,3,5]))
triangles = np.array(triangles)
for t in triangles:
points.append(((ratio+1)*points[t[0]]+ratio*points[t[1]])/(2*ratio+1.))
points.append(((ratio+1)*points[t[1]]+ratio*points[t[2]])/(2*ratio+1.))
points.append(((ratio+1)*points[t[2]]+ratio*points[t[0]])/(2*ratio+1.))
for t in triangles:
points.append(np.mean([points[p] for p in t],axis=1))
points = np.array(points)
points = center + radius*points
points = array_dict(dict(zip(np.arange(len(points)),points)))
faces = []
faces.append([6,16,27,19])
faces.append([9,7,18,22])
faces.append([12,10,21,25])
faces.append([15,13,24,28])
faces.append([8,11,14,17])
faces.append([20,23,26,29])
faces.append([6,16,8,11,7,18])
faces.append([9,7,11,14,10,21])
faces.append([12,10,14,17,13,24])
faces.append([15,13,17,8,16,27])
#faces.append([9,21,25,23,20,22])
faces.append([15,28,26,29,19,27])
faces.append([12,25,23,26,28,24])
faces.append([9,22,20,23,25,21])
faces.append([6,19,29,20,22,18])
# faces.append([9,7,11,14,10,21])
# faces.append([12,10,14,17,13,24])
# faces.append([15,13,17,8,16,27])
topomesh = poly_topomesh(faces,points)
return topomesh
[docs]def hexagonal_prism_topomesh(center=np.array([0,0,0],float),radius=10.,height=5.,angle=np.pi/6.):
points = {}
for p in range(6):
points[p] = [radius*np.cos(p*np.pi/3.+angle),radius*np.sin(p*np.pi/3.+angle),0]
top_hexagon_topomesh = poly_topomesh([[0,1,2,3,4,5]],points)
top_hexagon_topomesh.update_wisp_property('barycenter',0,center+top_hexagon_topomesh.wisp_property('barycenter',0).values()+np.array([0,0,height/2.]))
bottom_hexagon_topomesh = poly_topomesh([[0,1,2,3,4,5]],points)
bottom_hexagon_topomesh.update_wisp_property('barycenter',0,center+bottom_hexagon_topomesh.wisp_property('barycenter',0).values()-np.array([0,0,height/2.]))
top_positions = top_hexagon_topomesh.wisp_property('barycenter',0)
bottom_positions = bottom_hexagon_topomesh.wisp_property('barycenter',0)
side_topomeshes = []
for t_e, b_e in zip(top_hexagon_topomesh.wisps(1),bottom_hexagon_topomesh.wisps(1)):
top_vertices = list(top_hexagon_topomesh.borders(1,t_e))
bottom_vertices = list(bottom_hexagon_topomesh.borders(1,b_e))
points = {}
points[0] = top_positions[top_vertices[0]]
points[1] = top_positions[top_vertices[1]]
points[2] = bottom_positions[bottom_vertices[1]]
points[3] = bottom_positions[bottom_vertices[0]]
side_topomeshes.append(poly_topomesh([[0,1,2,3]], points))
topomesh,_ = append_topomesh(top_hexagon_topomesh,bottom_hexagon_topomesh)
for side_topomesh in side_topomeshes:
topomesh,_ = append_topomesh(topomesh,side_topomesh)
fuse_topomesh_identical_vertices(topomesh)
if not 0 in topomesh.wisps(3):
topomesh.add_wisp(3,0)
cells_to_remove = [c for c in topomesh.wisps(3) if c!=0]
for f in topomesh.wisps(2):
if not 0 in topomesh.regions(2,f):
topomesh.link(3,0,f)
for c in cells_to_remove:
topomesh.remove_wisp(3,c)
return topomesh
[docs]def sphere_topomesh(radius=1.0,center=np.zeros(3)):
from cellcomplex.utils import array_dict
from cellcomplex.property_topomesh.optimization import topomesh_triangle_split
ico = icosahedron_topomesh()
topomesh = topomesh_triangle_split(topomesh_triangle_split(ico))
positions = topomesh.wisp_property('barycenter',0)
new_positions = array_dict(center + radius*positions.values()/np.linalg.norm(positions.values(),axis=1)[:,np.newaxis],positions.keys())
topomesh.update_wisp_property('barycenter',0,new_positions)
return topomesh
[docs]def triangular_grid_topomesh(size=1, voxelsize=1., center=np.zeros(3)):
"""
"""
x = []
y = []
for n in range(size+1,2*size+2):
x += list((np.arange(n)-((n-1)/2.))*voxelsize*np.sqrt(3.))
y += [((2*size+1 - n))*voxelsize*3./2. for p in range(n)]
for n in range(size+1,2*size+1)[::-1]:
x += list((np.arange(n)-((n-1)/2.))*voxelsize*np.sqrt(3.))
y += [-((2*size+1 - n))*voxelsize*3./2. for p in range(n)]
x = np.array(x) + center[0]
y = np.array(y) + center[1]
z = np.zeros_like(x) + center[2]
positions = array_dict(dict(list(zip(list(range(len(x))),np.transpose([x,y,z])))))
triangles = positions.keys()[delaunay_triangulation(positions.values())]
topomesh = triangle_topomesh(triangles,positions)
return topomesh
[docs]def square_grid_topomesh(size=1, voxelsize=1., center=np.zeros(3)):
"""
"""
x = y = np.linspace(-size,size,2*size+1)*voxelsize
x += center[0]
y += center[1]
xx,yy = np.meshgrid(x,y)
zz = np.zeros_like(xx) + center[2]
square_bl = np.concatenate([np.arange(2*size+1)[:-1]+i*(2*size+1) for i in np.arange(2*size+1)[:-1]])
square_br = np.concatenate([np.arange(2*size+1)[1:]+i*(2*size+1) for i in np.arange(2*size+1)[:-1]])
square_tl = np.concatenate([np.arange(2*size+1)[:-1]+i*(2*size+1) for i in np.arange(2*size+1)[1:]])
square_tr = np.concatenate([np.arange(2*size+1)[1:]+i*(2*size+1) for i in np.arange(2*size+1)[1:]])
squares = np.transpose([square_bl,square_tl,square_tr,square_br])
positions = dict(list(zip(list(range(len(np.ravel(xx)))),np.transpose([np.ravel(xx),np.ravel(yy),np.ravel(zz)]))))
topomesh = quad_topomesh(squares,positions,faces_as_cells=True)
return topomesh
[docs]def hexagonal_grid_topomesh(size=1, voxelsize=1., center=np.zeros(3)):
"""
"""
x = []
y = []
p = []
for n in range(size+1,2*size+2):
for r in [0,1]:
p += list(len(p)+np.arange(n+r))
x += list((np.arange(n+r)-((n+r-1)/2.))*voxelsize*np.sqrt(3.))
row_y = (((2*size+1 - n)+1)//2)*3
row_side = (((2*size+1 - n)+1)%2)
y += list((row_y+(2*row_side-1)-(row_side==r)*(2*row_side-1)/2.)*np.ones(n+r)*voxelsize)
for n in range(size+1,2*size+2)[::-1]:
for r in [1,0]:
p += list(len(p)+np.arange(n+r))
x += list((np.arange(n+r)-((n+r-1)/2.))*voxelsize*np.sqrt(3.))
row_y = (((2*size+1 - n)+1)//2)*3
row_side = (((2*size+1 - n)+1)%2)
y += list(-(row_y+(2*row_side-1)-(row_side==r)*(2*row_side-1)/2.)*np.ones(n+r)*voxelsize)
x = np.array(x) + center[0]
y = np.array(y) + center[1]
z = np.zeros_like(x) + center[2]
positions = dict(list(zip(list(range(len(x))),np.transpose([x,y,z]))))
index_list = []
row_length_list = []
row_gaps = []
start = 0
for n in range(size+1,2*size+1):
index_list += list(np.arange(n)+start)
start += n+n+1
row_length_list += list(n*np.ones(n,int))
for h in range(n):
row_gaps += [[n+1,n+1,n+1,-(n+2),-(n+1)]]
for n in [2*size+1]:
index_list += list(np.arange(n)+start)
start += n+n+2
row_length_list += list(n*np.ones(n,int))
for h in range(n):
row_gaps += [[n+1,n+1,n,-(n+1),-(n+1)]]
for n in range(2*size,size,-1):
index_list += list(np.arange(n)+start)
start += n+n+3
row_length_list += list(n*np.ones(n,int))
for h in range(n):
row_gaps += [[n+2,n+1,n,-(n+1),-(n+1)]]
hexagons = []
for i,n,gaps in zip(index_list,row_length_list,row_gaps):
hexagon = [i]
for p, gap in zip(range(5),gaps):
hexagon += [hexagon[-1]+gap]
hexagons += [hexagon]
topomesh = poly_topomesh(hexagons,positions,faces_as_cells=True)
return topomesh
[docs]def circle_topomesh(radius = 1., n_points = 100, center=np.zeros(3,float)):
circle_thetas = np.linspace(-np.pi,np.pi-2*np.pi/float(n_points),n_points)
circle_points = center + np.transpose([radius*np.cos(circle_thetas),radius*np.sin(circle_thetas),np.zeros_like(circle_thetas)])
circle_positions = dict(zip(range(n_points),circle_points))
circle_vertices = np.arange(n_points)
return poly_topomesh([circle_vertices],circle_positions)
[docs]def circle_voronoi_topomesh(size = 1,voxelsize = 1.,circle_size = 100.,z_coef = 0.,iterations=None,return_triangulation=False):
n_cells = 3*size*(size-1)+1
radius = size*voxelsize
circle_thetas = np.linspace(-np.pi,np.pi-2*np.pi/float(circle_size),circle_size)
circle_points = np.transpose([radius*np.cos(circle_thetas),radius*np.sin(circle_thetas)])
cell_thetas = np.array([np.pi*np.random.randint(-180,180)/180. for c in range(n_cells)])
cell_distances = 0.5*radius*np.sqrt([np.random.rand() for c in range(n_cells)])
cell_points = np.transpose([cell_distances*np.cos(cell_thetas),cell_distances*np.sin(cell_thetas)])
omega_forces = dict(repulsion=0.5)
sigma_deformation = 2.*radius/float(n_cells)
if iterations is None:
iterations = n_cells//2
for iteration in range(iterations):
cell_to_cell_vectors = np.array([[p-q for q in cell_points] for p in cell_points])
cell_to_cell_distances = np.linalg.norm(cell_to_cell_vectors,axis=2)/radius
cell_to_circle_vectors = np.array([[p-q for q in circle_points] for p in cell_points])
cell_to_circle_distances = np.linalg.norm(cell_to_circle_vectors,axis=2)/radius
deformation_force = np.zeros_like(cell_points)
cell_repulsion_force = np.nansum(cell_to_cell_vectors/np.power(cell_to_cell_distances,3)[:,:,np.newaxis],axis=1)
circle_repulsion_force = np.nansum(cell_to_circle_vectors/np.power(cell_to_circle_distances,3)[:,:,np.newaxis],axis=1)
deformation_force += omega_forces['repulsion']*cell_repulsion_force
deformation_force += 1.5*(n_cells/float(circle_size))*omega_forces['repulsion']*circle_repulsion_force
deformation_force_amplitude = np.linalg.norm(deformation_force,axis=1)
deformation_force = np.minimum(1.0,sigma_deformation/deformation_force_amplitude)[:,np.newaxis] * deformation_force
cell_points += deformation_force
cell_points = np.minimum(1.0,radius/(np.linalg.norm(cell_points,axis=1)))[:,np.newaxis] * cell_points
all_positions = array_dict(np.transpose([np.concatenate([cell_points[:,0],circle_points[:,0]]),np.concatenate([cell_points[:,1],circle_points[:,1]]),np.zeros(int(np.round(n_cells+circle_size)))]),keys=np.arange(int(np.round(n_cells+circle_size))).astype(int))
triangles = all_positions.keys()[delaunay_triangulation(all_positions.values())]
radial_distances = np.linalg.norm(all_positions.values(),axis=1)
radial_z = z_coef*np.power(radial_distances/radius,2)
all_positions = array_dict(np.transpose([all_positions.values()[:,0],all_positions.values()[:,1],radial_z]),keys=all_positions.keys())
triangulation_topomesh = triangle_topomesh(triangles,all_positions)
cell_topomesh = dual_topomesh(triangulation_topomesh,2,vertex_positions='voronoi')
if return_triangulation:
return cell_topomesh, triangulation_topomesh
else:
return cell_topomesh
[docs]def vtk_ellipsoid_topomesh(ellipsoid_radius=50.0, ellipsoid_scales=[1,1,1], ellipsoid_axes=np.diag(np.ones(3)), ellipsoid_center=np.zeros(3), subdivisions=3):
"""
"""
import vtk
from cellcomplex.property_topomesh.utils.triangular_mesh_tools import vtk_polydata_to_triangular_mesh
ico = vtk.vtkPlatonicSolidSource()
ico.SetSolidTypeToIcosahedron()
ico.Update()
subdivide = vtk.vtkLoopSubdivisionFilter()
subdivide.SetNumberOfSubdivisions(subdivisions)
subdivide.SetInputConnection(ico.GetOutputPort())
subdivide.Update()
scale_transform = vtk.vtkTransform()
scale_factor = ellipsoid_radius/(np.sqrt(2)/2.)
scale_transform.Scale(scale_factor,scale_factor,scale_factor)
ellipsoid_sphere = vtk.vtkTransformPolyDataFilter()
ellipsoid_sphere.SetInputData(subdivide.GetOutput())
ellipsoid_sphere.SetTransform(scale_transform)
ellipsoid_sphere.Update()
ellipsoid_transform = vtk.vtkTransform()
axes_transform = vtk.vtkLandmarkTransform()
source_points = vtk.vtkPoints()
source_points.InsertNextPoint([1,0,0])
source_points.InsertNextPoint([0,1,0])
source_points.InsertNextPoint([0,0,1])
target_points = vtk.vtkPoints()
target_points.InsertNextPoint(ellipsoid_axes[0])
target_points.InsertNextPoint(ellipsoid_axes[1])
target_points.InsertNextPoint(ellipsoid_axes[2])
axes_transform.SetSourceLandmarks(source_points)
axes_transform.SetTargetLandmarks(target_points)
axes_transform.SetModeToRigidBody()
axes_transform.Update()
ellipsoid_transform.SetMatrix(axes_transform.GetMatrix())
ellipsoid_transform.Scale(ellipsoid_scales[0],ellipsoid_scales[1],ellipsoid_scales[2])
center_transform = vtk.vtkTransform()
center_transform.Translate(ellipsoid_center[0],ellipsoid_center[1],ellipsoid_center[2])
center_transform.Concatenate(ellipsoid_transform)
ellipsoid_ellipsoid = vtk.vtkTransformPolyDataFilter()
ellipsoid_ellipsoid.SetInputData(ellipsoid_sphere.GetOutput())
ellipsoid_ellipsoid.SetTransform(center_transform)
ellipsoid_ellipsoid.Update()
ellipsoid_mesh = vtk_polydata_to_triangular_mesh(ellipsoid_ellipsoid.GetOutput())
triangles = np.array([t for t in ellipsoid_mesh.triangles.values()])
topomesh = triangle_topomesh(triangles,ellipsoid_mesh.points)
return topomesh
[docs]def hexagonal_prism_layered_grid_topomesh(size=1, n_layers=1,center=np.array([0,0,0]), radius=10, height_factor=2/3., triangulate=True, smoothing=False):
topomesh = hexagonal_prism_topomesh(center=center,radius=radius,height=height_factor*radius,angle=np.pi/6.)
for layer in range(n_layers):
for row in np.arange(-size,size+1):
for col in (np.arange(2*size+1 - np.abs(row) - layer) - (size-layer/2.-(np.abs(row))/2.)):
layer_row = row-layer/3.
if ((layer_row!=0) or (col!=0) or (layer!=0))&(np.abs(layer_row)<=(size-(layer-1)/2.)):
topomesh,_ = append_topomesh(topomesh, hexagonal_prism_topomesh(center=np.array([2.*col*radius*np.sqrt(3.)/2.,2.*layer_row*radius*3./4.,-layer*height_factor*radius]),radius=radius,height=height_factor*radius,angle=np.pi/6.))
for layer in range(n_layers):
for row in np.arange(-size,size+1):
for col in (np.arange(2*size+1 - np.abs(row) - layer) - (size-layer/2.-(np.abs(row))/2.)):
layer_row = row-layer/3.
print(layer,layer_row, col)
if triangulate:
topomesh = star_interface_topomesh(topomesh)
fuse_topomesh_identical_vertices(topomesh)
if triangulate and smoothing:
topomesh = topomesh_triangle_split(topomesh)
for cid in topomesh.wisps(3):
for n_cid in topomesh.border_neighbors(3,cid):
if (n_cid<cid) and (not (n_cid,cid) in topomesh._interface[3].values()):
iid = topomesh._interface[3].add((n_cid,cid),None)
compute_topomesh_property(topomesh,'vertices',degree=1)
compute_topomesh_property(topomesh,'vertices',degree=2)
compute_topomesh_property(topomesh,'vertices',degree=3)
if triangulate and smoothing:
compute_topomesh_triangle_properties(topomesh)
property_topomesh_vertices_deformation(topomesh,omega_forces=dict(taubin_smoothing=0.65,planarization=1),iterations=100)
# property_topomesh_vertices_deformation(topomesh,omega_forces=dict(taubin_smoothing=0.65,regularization=0.01,planarization=1),iterations=50)
# topomesh = property_topomesh_isotropic_remeshing(topomesh,maximal_length=3.,iterations=10)
return topomesh
[docs]def truncated_octahedra_layered_grid_topomesh(size=1, n_layers=1, center=np.array([0,0,0]), radius=10, triangulate=True, smoothing=False):
topomesh = truncated_octahedron_topomesh(center=center,radius=radius)
for layer in range(n_layers):
for row in np.arange(-size,size+1):
for col in (np.arange(2*size+1 - np.abs(row) - layer) - (size-layer/2.-(np.abs(row))/2.)):
print(layer,row, col)
if (row!=0) or (col!=0) or (layer!=0):
topomesh,_ = append_topomesh(topomesh, truncated_octahedron_topomesh(center=np.array([4.*col*radius*np.sqrt(2)/3,4.*row*radius/3,-4.*layer*radius/3]),radius=radius))
fuse_topomesh_identical_vertices(topomesh)
if triangulate:
topomesh = star_interface_topomesh(topomesh)
if triangulate and smoothing:
topomesh = topomesh_triangle_split(topomesh)
for cid in topomesh.wisps(3):
for n_cid in topomesh.border_neighbors(3,cid):
if (n_cid<cid) and (not (n_cid,cid) in topomesh._interface[3].values()):
iid = topomesh._interface[3].add((n_cid,cid),None)
compute_topomesh_property(topomesh,'vertices',degree=1)
compute_topomesh_property(topomesh,'vertices',degree=2)
compute_topomesh_property(topomesh,'vertices',degree=3)
if triangulate and smoothing:
compute_topomesh_triangle_properties(topomesh)
property_topomesh_vertices_deformation(topomesh,omega_forces=dict(taubin_smoothing=0.65,planarization=1),iterations=100)
# property_topomesh_vertices_deformation(topomesh,omega_forces=dict(taubin_smoothing=0.65,regularization=0.01,planarization=1),iterations=50)
# topomesh = property_topomesh_isotropic_remeshing(topomesh,maximal_length=3.,iterations=10)
return topomesh
