# -*- 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/
#
###############################################################################
"""
=========================
PropertyTopomesh Analysis
=========================
Implementing functions to compute geometrical and topological properties on
a PropertyTopomesh, generally assuming it represents a triangular mesh (but
not necessarily a 2-manifold one). Powerful computation of the values of the
property for all the mesh elements at once. Convenient functions such as edge
lengths, cell volumes, surface curvatures, vertex normals, face areas, and
so on.
.. currentmodule:: cellcomplex.property_topomesh.analysis
.. autosummary::
:toctree: generated/
compute_topomesh_property
compute_topomesh_triangle_properties
compute_topomesh_vertex_property_from_faces
compute_topomesh_cell_property_from_faces
topomesh_property_gaussian_filtering
"""
from copy import copy, deepcopy
from time import time
from array import array
import os
import sys
import pickle
import numpy as np
from scipy import ndimage as nd
from scipy.cluster.vq import kmeans, vq
from cellcomplex.property_topomesh import PropertyTopomesh
from cellcomplex.utils import array_dict
[docs]def compute_topomesh_property(topomesh, property_name, degree=0, positions=None, normal_method="density", object_positions=None, object_radius=10., verbose=False):
"""Compute a property over the elements of a PropertyTopomesh.
The function computes and fills a property of a PropertyTopomesh passed as
argument. The given property is computed for all elements of the specified
degree and stored as a dictionary in the PropertyTopomesh structure.
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to compute the property.
property_name : str
The name of the property to compute, among the following ones:
* *barycenter* (degree : [0, 1, 2, 3])
* *vertices* (degree : [0, 1, 2, 3])
* *triangles* (degree : [0, 1, 2, 3])
* *cells* (degree : [0, 1, 2, 3])
* *length* (degree : [1])
* *area* (degree : [2])
* *volume* (degree : [3])
* *normal* (degree : [0, 2])
* *skewness* (degree : [2])
* *roughness* (degree : [0])
* *vertex_neighbors* (degree : [2])
degree : int
The degree of the elements on which to compute the property.
positions : dict
A position dictionary to use if the property ('barycenter',0) is empty.
normal_method : str
The method used to re-orient the face normals (default is *density*)
* *barycenter*: oriented using the direction of the mesh center
* *density*: oriented using the gradient of object density
* *orientation*: consistently oriented using topological propagation
object_positions : dict
The position of the object(s) represented by the mesh.
* used only for property ('normal',2)
object_radius : float
The radius of the object(s) represented by the mesh.
* used only for property ('normal',2)
verbose : bool
Whether to display or not information on computed properties.
Returns
--------
None
Note
-------
The PropertyTopomesh passed as argument is updated.
Example
-------
>>> from cellcomplex.property_topomesh.example_topomesh import square_topomesh
>>> from cellcomplex.property_topomesh.analysis import compute_topomesh_property
>>> topomesh = square_topomesh(side_length=1)
>>> compute_topomesh_property(topomesh,'length',1)
>>> print topomesh.wisp_property('length',1)
{0: 1.0, 1: 1.0, 2: 1.41421356237, 3: 1.0, 4: 1.0}
"""
if positions is None:
positions = topomesh.wisp_property('barycenter',degree=0)
start_time = time()
if verbose:
print("--> Computing ",property_name," property (",degree,")")
if property_name == 'barycenter':
if not isinstance(positions,array_dict):
positions = array_dict(positions)
if not 'barycenter' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('barycenter',degree=degree)
if degree == 0:
topomesh.update_wisp_property('barycenter',degree=degree,values=positions,keys=np.array(list(topomesh.wisps(degree))))
elif degree == 1:
compute_topomesh_property(topomesh,'vertices',degree=degree)
element_vertices = topomesh.wisp_property('vertices',degree).values(np.array(list(topomesh.wisps(degree))))
element_vertices_number = np.array(list(map(len,element_vertices)))
barycentrable_elements = np.array(list(topomesh.wisps(degree)))[np.where(element_vertices_number>0)]
element_vertices = topomesh.wisp_property('vertices',degree).values(barycentrable_elements)
vertices_positions = positions.values(element_vertices)
if vertices_positions.dtype == np.dtype('O'):
barycenter_positions = np.array([np.mean(p,axis=0) for p in vertices_positions])
else:
barycenter_positions = np.mean(vertices_positions,axis=1)
topomesh.update_wisp_property('barycenter',degree=degree,values=barycenter_positions,keys=barycentrable_elements)
elif degree == 2:
if not topomesh.has_wisp_property('barycenter',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=1)
if not topomesh.has_wisp_property('length',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'length',degree=1)
compute_topomesh_property(topomesh,'edges',degree=degree)
element_edges = topomesh.wisp_property('edges',degree).values(np.array(list(topomesh.wisps(degree))))
element_edges_number = np.array(list(map(len,element_edges)))
barycentrable_elements = np.array(list(topomesh.wisps(degree)))[np.where(element_edges_number>0)]
element_edges = topomesh.wisp_property('edges',degree).values(barycentrable_elements)
edge_positions = topomesh.wisp_property('barycenter',1).values(element_edges)
edge_weights = topomesh.wisp_property('length',1).values(element_edges)
if edge_positions.dtype == np.dtype('O'):
barycenter_positions = np.array([np.sum(w[:,np.newaxis]*p,axis=0)/np.sum(w) for p,w in zip(edge_positions,edge_weights)])
else:
barycenter_positions = np.sum(edge_weights[:,:,np.newaxis]*edge_positions,axis=1)/np.sum(edge_weights,axis=1)[:,np.newaxis]
topomesh.update_wisp_property('barycenter',degree=degree,values=barycenter_positions,keys=barycentrable_elements)
elif degree == 3:
if not topomesh.has_wisp_property('barycenter', degree=2, is_computed=True):
compute_topomesh_property(topomesh, 'barycenter', degree=2)
if not topomesh.has_wisp_property('area', degree=2, is_computed=True):
compute_topomesh_property(topomesh, 'area', degree=2)
compute_topomesh_property(topomesh, 'faces', degree=degree)
element_faces = topomesh.wisp_property('faces', degree).values(np.array(list(topomesh.wisps(degree))))
element_faces_number = np.array(list(map(len, element_faces)))
barycentrable_elements = np.array(list(topomesh.wisps(degree)))[np.where(element_faces_number > 0)]
element_faces = topomesh.wisp_property('faces', degree).values(barycentrable_elements)
face_positions = topomesh.wisp_property('barycenter', 2).values(element_faces)
face_weights = topomesh.wisp_property('area', 2).values(element_faces)
if face_positions.dtype == np.dtype('O'):
barycenter_positions = np.array([np.sum(w[:, np.newaxis] * p, axis=0) / np.sum(w) for p, w in zip(face_positions, face_weights)])
else:
barycenter_positions = np.sum(face_weights[:, :, np.newaxis] * face_positions, axis=1) / np.sum(face_weights, axis=1)[:, np.newaxis]
topomesh.update_wisp_property('barycenter', degree=degree, values=barycenter_positions, keys=barycentrable_elements)
if property_name == 'vertices':
if not 'vertices' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('vertices',degree=degree)
if degree == 0:
topomesh.update_wisp_property('vertices',degree=degree,values=np.array(list(topomesh.wisps(degree))).astype(int),keys=np.array(list(topomesh.wisps(degree))))
else:
if verbose: print(list(topomesh.wisps(degree)))
topomesh.update_wisp_property('vertices',degree=degree,values=np.array([np.unique(list(topomesh.borders(degree,w,degree))).astype(int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
if degree == 1:
topomesh.update_wisp_property('borders',degree=degree,values=np.array([np.unique(list(topomesh.borders(degree,w,degree))).astype(int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'edges':
if not 'edges' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('edges',degree=degree)
if degree < 1:
topomesh.update_wisp_property('edges',degree=degree,values=np.array([np.array(list(topomesh.regions(degree,w,1-degree))).astype(int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
elif degree > 1:
topomesh.update_wisp_property('edges',degree=degree,values=np.array([np.array(list(topomesh.borders(degree,w,degree-1))).astype(int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
else:
topomesh.update_wisp_property('edges',degree=degree,values=np.array(list(topomesh.wisps(degree))).astype(int),keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'faces':
if not 'faces' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('faces',degree=degree)
if degree < 2:
topomesh.update_wisp_property('faces',degree=degree,values=np.array([np.array(list(topomesh.regions(degree,w,2-degree))).astype(int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
elif degree > 2:
topomesh.update_wisp_property('faces',degree=degree,values=np.array([np.array(list(topomesh.borders(degree,w,degree-2))).astype(int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
else:
topomesh.update_wisp_property('faces',degree=degree,values=np.array(list(topomesh.wisps(degree))).astype(int),keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'cells':
if not 'cells' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('cells',degree=degree)
if degree == topomesh.degree():
topomesh.update_wisp_property('cells',degree=degree,values=np.array(list(topomesh.wisps(degree))).astype(int),keys=np.array(list(topomesh.wisps(degree))))
else:
topomesh.update_wisp_property('cells',degree=degree,values=np.array([[int(c) for c in topomesh.regions(degree,w,topomesh.degree()-degree)] for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'borders':
assert degree>0
if not 'borders' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('borders',degree=degree)
topomesh.update_wisp_property('borders',degree=degree,values=np.array([np.array(list(topomesh.borders(degree,w)),int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'neighbors':
if not 'neighbors' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('neighbors',degree=degree)
if degree>0:
topomesh.update_wisp_property('neighbors',degree=degree,values=np.array([np.unique(list(topomesh.border_neighbors(degree,w))).astype(int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
else:
topomesh.update_wisp_property('neighbors',degree=degree,values=np.array([np.unique(list(topomesh.region_neighbors(degree,w))).astype(int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'regions':
assert degree<topomesh.degree()
if not 'regions' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('regions',degree=degree)
topomesh.update_wisp_property('regions',degree=degree,values=np.array([np.array(list(topomesh.regions(degree,w))).astype(int) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'valency':
assert degree<topomesh.degree()
if not 'valency' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('valency',degree=degree)
# if not topomesh.has_wisp_property('neighbors',degree=degree,is_computed=True):
# compute_topomesh_property(topomesh,'neighbors',degree=degree)
# topomesh.update_wisp_property('valency',degree=degree,values=np.array([len(topomesh.wisp_property('neighbors',degree)[w]) for w in topomesh.wisps(degree)]),keys=np.array(list(topomesh.wisps(degree))))
topomesh.update_wisp_property('valency',degree=degree,values=np.array([len(list(topomesh.region_neighbors(degree,w))) for w in topomesh.wisps(degree)]).astype(int),keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'oriented_borders':
if degree == 2:
oriented_borders = []
oriented_border_orientations = []
for t,fid in enumerate(topomesh.wisps(2)):
face_eids = np.array(list(topomesh.borders(2,fid)))
#face_edges = np.array([list(topomesh.borders(1,eid)) for eid in face_eids])
oriented_face_eids = [face_eids[0]]
oriented_face_eid_orientations = [1]
while len(oriented_face_eids) < len(face_eids):
current_eid = oriented_face_eids[-1]
current_eid_orientation = oriented_face_eid_orientations[-1]
if current_eid_orientation == 1:
start_pid, end_pid = topomesh.borders(1,current_eid)
else:
end_pid, start_pid = topomesh.borders(1,current_eid)
candidate_eids = set(list(topomesh.regions(0,end_pid))).intersection(set(list(face_eids))).difference({current_eid})
if len(oriented_face_eids)>0:
candidate_eids = candidate_eids.difference(set(oriented_face_eids))
if len(candidate_eids)>0:
next_eid = list(candidate_eids)[0]
else:
next_eid = list(set(list(face_eids)).difference(set(oriented_face_eids)))[0]
oriented_face_eids += [next_eid]
if end_pid == list(topomesh.borders(1,next_eid))[0]:
oriented_face_eid_orientations += [1]
else:
oriented_face_eid_orientations += [-1]
# print oriented_face_eids," (",oriented_face_eid_orientations,")"
oriented_borders += [np.array(oriented_face_eids)]
oriented_border_orientations += [np.array(oriented_face_eid_orientations)]
# raw_input()
topomesh.update_wisp_property('oriented_borders',degree=2,values=dict(zip(topomesh.wisps(2),[[b,o] for b,o in zip(oriented_borders,oriented_border_orientations)])))
elif degree == 3:
compute_topomesh_property(topomesh,'oriented_borders',2)
compute_topomesh_property(topomesh,'oriented_vertices',2)
compute_topomesh_property(topomesh,'barycenter',2)
compute_topomesh_property(topomesh,'barycenter',3)
triangle_edge_list = np.array([[0, 1],[1, 2],[2, 0]])
oriented_borders = []
oriented_border_orientations = []
oriented_border_components = []
oriented_face_normals = []
for f,fid in enumerate(topomesh.wisps(2)):
face_vertices = topomesh.wisp_property('oriented_vertices',2)[fid]
if len(face_vertices) == 3:
vertices_positions = topomesh.wisp_property('barycenter',0).values(face_vertices)
face_normal = np.cross(vertices_positions[1]-vertices_positions[0],vertices_positions[2]-vertices_positions[0])
else:
face_edges = np.transpose(np.array([face_vertices,list(face_vertices[1:])+[face_vertices[0]]]))
face_barycenter = topomesh.wisp_property('barycenter',2)[fid]
vertices_positions = np.array([np.concatenate([p,[face_barycenter]]) for p in topomesh.wisp_property('barycenter',0).values(face_edges)])
face_triangle_normals = np.cross(vertices_positions[:,1]-vertices_positions[:,0],vertices_positions[:,2]-vertices_positions[:,0])
face_triangle_normals = face_triangle_normals/np.linalg.norm(face_triangle_normals,axis=1)[:,np.newaxis]
face_triangle_lengths = np.linalg.norm(vertices_positions[...,triangle_edge_list,1] - vertices_positions[...,triangle_edge_list,0],axis=2)
face_triangle_perimeters = face_triangle_lengths.sum(axis=1)
face_triangle_areas = np.sqrt(np.maximum(0,(face_triangle_perimeters/2.0)*(face_triangle_perimeters/2.0-face_triangle_lengths[...,0])*(face_triangle_perimeters/2.0-face_triangle_lengths[...,1])*(face_triangle_perimeters/2.0-face_triangle_lengths[...,2])))
face_normal = (face_triangle_normals*face_triangle_areas[:,np.newaxis]).sum(axis=0)/face_triangle_areas.sum()
face_normal = face_normal/np.linalg.norm(face_normal)
oriented_face_normals += [face_normal]
topomesh.update_wisp_property('normal',degree=2,values=np.array(oriented_face_normals),keys=np.array(list(topomesh.wisps(2))))
cell_face_component = 0
for c,cid in enumerate(topomesh.wisps(3)):
cell_fids = np.array(list(topomesh.borders(3,cid))).astype(int)
# print cell_fids
cell_fid_normals = topomesh.wisp_property('normal',2).values(cell_fids)
cell_fid_barycenters = topomesh.wisp_property('barycenter',2).values(cell_fids)
cell_barycenter = topomesh.wisp_property('barycenter',3)[cid]
cell_fid_dot_product = np.einsum('...ij,...ij->...i',cell_fid_normals,cell_fid_barycenters-cell_barycenter[np.newaxis])
cell_fid_dot_product_signs = np.array(np.sign(cell_fid_dot_product),int)
cell_fid_dot_product_signs[cell_fid_dot_product_signs==0] = 2*(cell_fid_normals[cell_fid_dot_product_signs==0][:,2]>0) - 1
cell_fid_orientation = array_dict(cell_fid_dot_product_signs,cell_fids)
cell_fid_dot_product = array_dict(cell_fid_dot_product,cell_fids)
start_fid = cell_fids[np.argmax(np.abs(cell_fid_dot_product.values(cell_fids)))]
#start_fid = cell_fids[0]
oriented_cell_fid_orientations = {}
oriented_cell_fid_orientations[start_fid] = cell_fid_orientation[start_fid]
cell_face_component = cell_face_component+1
fid_components = {}
fid_components[start_fid] = cell_face_component
next_fid_queue = []
current_fid_queue = []
transition_eid_queue = []
eid_orientations = array_dict(*topomesh.wisp_property('oriented_borders',2)[start_fid][::-1])
for eid in topomesh.wisp_property('oriented_borders',2)[start_fid][0]:
candidate_fids = set(list(topomesh.regions(1,eid))).intersection(set(list(cell_fids))).difference(set(oriented_cell_fid_orientations.keys()))
if len(candidate_fids)>0:
for fid in candidate_fids:
next_fid_queue.insert(0,fid)
current_fid_queue.insert(0,start_fid)
transition_eid_queue.insert(0,eid)
# print oriented_cell_fid_orientations
# print next_fid_queue
# print len(oriented_cell_fid_orientations.keys())," / ",len(cell_fids)," [",len(next_fid_queue),"]"
# raw_input()
while len(oriented_cell_fid_orientations.keys()) < len(cell_fids):
while len(next_fid_queue) > 0:
next_fid = next_fid_queue.pop()
current_fid = current_fid_queue.pop()
transition_eid = transition_eid_queue.pop()
current_face_eid_orientations = array_dict(*topomesh.wisp_property('oriented_borders',2)[current_fid][::-1])
current_fid_orientation = oriented_cell_fid_orientations[current_fid]
next_face_eid_orientations = array_dict(*topomesh.wisp_property('oriented_borders',2)[next_fid][::-1])
# print list(topomesh.borders(2,current_fid))
# print current_face_eid_orientations
# print list(topomesh.borders(2,next_fid))
# print next_face_eid_orientations
next_fid_orientation = - current_face_eid_orientations[transition_eid]*current_fid_orientation*next_face_eid_orientations[transition_eid]
oriented_cell_fid_orientations[next_fid] = next_fid_orientation
fid_components[next_fid] = cell_face_component
# print list(topomesh.border_neighbors(2,next_fid))
for eid in topomesh.wisp_property('oriented_borders',2)[next_fid][0]:
candidate_fids = set(list(topomesh.regions(1,eid))).intersection(set(list(cell_fids)))
candidate_fids = candidate_fids.difference(set(oriented_cell_fid_orientations.keys()))
candidate_fids = candidate_fids.difference(set(next_fid_queue))
if len(candidate_fids)>0:
for fid in candidate_fids:
next_fid_queue.insert(0,fid)
current_fid_queue.insert(0,next_fid)
transition_eid_queue.insert(0,eid)
# print oriented_cell_fid_orientations
# next_fid_queue
# print "Cell ",cid," : ", len(oriented_cell_fid_orientations.keys())," / ",len(cell_fids)," [",len(next_fid_queue),"]"
# raw_input()
if len(oriented_cell_fid_orientations.keys()) < len(cell_fids):
remaining_cell_fids = np.array(list(set(cell_fids).difference(set(oriented_cell_fid_orientations.keys()))))
start_fid = remaining_cell_fids[np.argmax(np.abs(cell_fid_dot_product.values(remaining_cell_fids)))]
oriented_cell_fid_orientations[start_fid] = cell_fid_orientation[start_fid]
cell_face_component = cell_face_component+1
fid_components[start_fid] = cell_face_component
eid_orientations = array_dict(*topomesh.wisp_property('oriented_borders',2)[start_fid][::-1])
for eid in topomesh.wisp_property('oriented_borders',2)[start_fid][0]:
candidate_fids = set(list(topomesh.regions(1,eid))).intersection(set(list(cell_fids))).difference(set(oriented_cell_fid_orientations.keys()))
if len(candidate_fids)>0:
for fid in candidate_fids:
next_fid_queue.insert(0,fid)
current_fid_queue.insert(0,start_fid)
transition_eid_queue.insert(0,eid)
# print oriented_cell_fid_orientations
# print next_fid_queue
# print len(oriented_cell_fid_orientations.keys())," / ",len(cell_fids)," [",len(next_fid_queue),"]"
# raw_input()
# oriented_cell_fids = [cell_fids[0]]
# oriented_cell_fid_orientations = [cell_fid_orientation[cell_fids[0]]]
# start_eid = topomesh.wisp_property('oriented_borders',2)[cell_fids[0]][0][0]
# face_gaps = 0
# while len(oriented_cell_fids) < len(cell_fids):
# current_fid = oriented_cell_fids[-1]
# current_fid_orientation = oriented_cell_fid_orientations[-1]
# face_eids = np.array(topomesh.wisp_property('oriented_borders',2)[current_fid][0])
# face_eid_orientation = array_dict(np.array(topomesh.wisp_property('oriented_borders',2)[current_fid][1]),face_eids)
# candidate_fids = set()
# current_eid = start_eid
# n_tries = 0
# while (len(list(candidate_fids))==0) and n_tries<len(face_eids):
# candidate_fids = set(list(topomesh.regions(1,current_eid))).intersection(set(list(cell_fids))).difference(set(oriented_cell_fids))
# # print current_eid,"(",face_eids,") -> ",candidate_fids
# previous_eid = current_eid
# current_eid = face_eids[(np.where(face_eids == current_eid)[0][0] + current_fid_orientation) % len(face_eids)]
# n_tries += 1
# if n_tries == len(face_eids):
# index = 0
# remaining_triangles = len(list(set(cell_fids).difference(oriented_cell_fids)))
# next_fid = list(set(cell_fids).difference(oriented_cell_fids))[index]
# edge_index = 0
# start_eid = topomesh.wisp_property('oriented_borders',2)[next_fid][0][edge_index]
# opposite_fid = list(set(topomesh.regions(1,start_eid)).difference({next_fid}))[0] if topomesh.nb_regions(1,start_eid)>1 else -1
# while (not opposite_fid in oriented_cell_fids) and (edge_index < 2):
# edge_index = edge_index+1
# start_eid = topomesh.wisp_property('oriented_borders',2)[next_fid][0][edge_index]
# opposite_fid = list(set(topomesh.regions(1,start_eid)).difference({next_fid}))[0] if topomesh.nb_regions(1,start_eid)>1 else -1
# while (not opposite_fid in oriented_cell_fids) and (index < remaining_triangles):
# next_fid = list(set(cell_fids).difference(oriented_cell_fids))[index]
# index = index+1
# edge_index = 0
# start_eid = topomesh.wisp_property('oriented_borders',2)[next_fid][0][edge_index]
# opposite_fid = list(set(topomesh.regions(1,start_eid)).difference({next_fid}))[0] if topomesh.nb_regions(1,start_eid)>1 else -1
# while (not opposite_fid in oriented_cell_fids) and (edge_index < 2):
# edge_index = edge_index+1
# start_eid = topomesh.wisp_property('oriented_borders',2)[next_fid][0][edge_index]
# opposite_fid = list(set(topomesh.regions(1,start_eid)).difference({next_fid}))[0] if topomesh.nb_regions(1,start_eid)>1 else -1
# if index < len(list(set(cell_fids).difference(oriented_cell_fids))):
# #next_fid_orientation = cell_fid_orientation[next_fid]
# #next_fid_orientation = oriented_cell_fid_orientations
# opposite_face_eids = np.array(topomesh.wisp_property('oriented_borders',2)[opposite_fid][0])
# opposite_face_eid_orientation = array_dict(np.array(topomesh.wisp_property('oriented_borders',2)[opposite_fid][1]),opposite_face_eids)
# next_face_eids = np.array(topomesh.wisp_property('oriented_borders',2)[next_fid][0])
# next_face_eid_orientation = array_dict(np.array(topomesh.wisp_property('oriented_borders',2)[next_fid][1]),next_face_eids)
# opposite_fid_orientation = np.array(oriented_cell_fid_orientations)[np.where(oriented_cell_fids==opposite_fid)]
# next_fid_orientation = - opposite_face_eid_orientation[start_eid] * current_fid_orientation * next_face_eid_orientation[start_eid]
# else:
# index = 0
# next_fid = list(set(cell_fids).difference(oriented_cell_fids))[index]
# start_eid = topomesh.wisp_property('oriented_borders',2)[next_fid][0][0]
# next_fid_orientation = cell_fid_orientation[next_fid]
# face_gaps += 1
# else:
# next_fid = int(list(candidate_fids)[0])
# start_eid = previous_eid
# next_face_eids = np.array(topomesh.wisp_property('oriented_borders',2)[next_fid][0])
# next_face_eid_orientation = array_dict(np.array(topomesh.wisp_property('oriented_borders',2)[next_fid][1]),next_face_eids)
# next_fid_orientation = - face_eid_orientation[start_eid]*current_fid_orientation * next_face_eid_orientation[start_eid]
# oriented_cell_fids += [next_fid]
# oriented_cell_fid_orientations += [next_fid_orientation]
# # print cell_fids," -> ",oriented_cell_fids,"(",face_gaps,"face gaps)"
# # print cell_fid_orientation.values(oriented_cell_fids) - np.array(oriented_cell_fid_orientations)
oriented_borders += [oriented_cell_fid_orientations.keys()]
oriented_face_normals = topomesh.wisp_property('normal',2).values(list(oriented_cell_fid_orientations.keys()))
oriented_face_oriented_normals = np.array([oriented_cell_fid_orientations[b] for b in oriented_cell_fid_orientations.keys()])[:,np.newaxis]*oriented_face_normals
oriented_face_barycenters = topomesh.wisp_property('barycenter',2).values(list(oriented_cell_fid_orientations.keys()))
oriented_face_dot_products = np.einsum('...ij,...ij->...i',oriented_face_oriented_normals,oriented_face_barycenters-cell_barycenter[np.newaxis])
oriented_face_dot_product_signs = np.array(np.sign(oriented_face_dot_products),int)
# oriented_border_orientations += [[int(np.median(oriented_face_dot_product_signs)*oriented_cell_fid_orientations[b]) for b in oriented_cell_fid_orientations.keys()]]
# print(oriented_face_dot_product_signs)
# print(oriented_cell_fid_orientations)
#
median_sign = np.median(oriented_face_dot_product_signs)
median_sign += int(median_sign==0)
oriented_border_orientations += [[int(median_sign*oriented_cell_fid_orientations[b]) for b in oriented_cell_fid_orientations.keys()]]
oriented_border_components += [[fid_components[b] for b in oriented_cell_fid_orientations.keys()]]
# oriented_border_orientations += [[oriented_cell_fid_orientations[b] for b in oriented_cell_fid_orientations.keys()]]
topomesh.update_wisp_property('oriented_borders',degree=3,values=np.array([(b,o) for b,o in zip(oriented_borders,oriented_border_orientations)]),keys=np.array(list(topomesh.wisps(3))))
topomesh.update_wisp_property('oriented_border_components',degree=3,values=np.array([(b,c) for b,c in zip(oriented_borders,oriented_border_components)]),keys=np.array(list(topomesh.wisps(3))))
# oriented_face_pids = [list(topomesh.borders(1,eid))[0 if ori==1 else 1] for eid,ori in zip(oriented_face_eids,oriented_face_eid_orientations)]
# print oriented_face_pids
# raw_input()
# face_pids = np.array(list(topomesh.borders(2,fid,2)))
# face_edges = np.array([list(topomesh.borders(1,eid)) for eid in topomesh.borders(2,fid)])
# oriented_face_pids = [face_pids[0]]
# while len(oriented_face_pids) < len(face_pids):
# current_pid = oriented_face_pids[-1]
# pid_edges = face_edges[np.where(face_edges==current_pid)[0]]
# candidate_pids = set(list(np.unique(pid_edges))).difference({current_pid})
# if len(oriented_face_pids)>1:
# candidate_pids = candidate_pids.difference({oriented_face_pids[-2]})
# oriented_face_pids += [list(candidate_pids)[0]]
# print oriented_face_pids
# raw_input()
if property_name == 'oriented_vertices':
if degree == 2:
compute_topomesh_property(topomesh,'oriented_borders',2)
oriented_vertices = []
for t,fid in enumerate(topomesh.wisps(2)):
oriented_face_eids = topomesh.wisp_property('oriented_borders',2)[fid][0]
oriented_face_eid_orientations = topomesh.wisp_property('oriented_borders',2)[fid][1]
oriented_face_pids = [list(topomesh.borders(1,eid))[0 if ori==1 else 1] for eid,ori in zip(oriented_face_eids,oriented_face_eid_orientations)]
oriented_vertices += [oriented_face_pids]
topomesh.update_wisp_property('oriented_vertices',degree=2,values=np.array(oriented_vertices),keys=np.array(list(topomesh.wisps(2))))
# face_pids = np.array(list(topomesh.borders(2,fid,2)))
# face_edges = np.array([list(topomesh.borders(1,eid)) for eid in topomesh.borders(2,fid)])
# oriented_face_pids = [face_pids[0]]
# while len(oriented_face_pids) < len(face_pids):
# current_pid = oriented_face_pids[-1]
# pid_edges = face_edges[np.where(face_edges==current_pid)[0]]
# candidate_pids = set(list(np.unique(pid_edges))).difference({current_pid})
# if len(oriented_face_pids)>1:
# candidate_pids = candidate_pids.difference({oriented_face_pids[-2]})
# oriented_face_pids += [list(candidate_pids)[0]]
# if property_name == 'curvature':
# assert degree==0
# if not 'curvature' in topomesh.wisp_property_names(degree):
# topomesh.add_wisp_property('curvature',degree=degree)
# if not topomesh.has_wisp_property('vertices',degree=1,is_computed=True):
# compute_topomesh_property(topomesh,'vertices',degree=1)
# if not topomesh.has_wisp_property('vertices',degree=3,is_computed=True):
# compute_topomesh_property(topomesh,'vertices',degree=3)
# if not topomesh.has_wisp_property('edges',degree=3,is_computed=True):
# compute_topomesh_property(topomesh,'edges',degree=3)
# if not topomesh.has_wisp_property('cells',degree=0,is_computed=True):
# compute_topomesh_property(topomesh,'cells',degree=0)
# if not topomesh.has_wisp_property('neighbors',degree=0,is_computed=True):
# compute_topomesh_property(topomesh,'neighbors',degree=0)
# vertex_cells = topomesh.wisp_property('cells',0).values()
# vertex_cell_neighbors = np.array([[np.intersect1d(topomesh.wisp_property('vertices',3)[c],topomesh.wisp_property('neighbors',0)[v]) for c in topomesh.wisp_property('cells',0)[v]] for v in topomesh.wisps(0)])
# vertex_cell_valency = np.array([[len(n) for n in vertex_cell_neighbors[v]] for v in topomesh.wisps(0)])
# vertex_cell_edge_vectors = np.array([[topomesh.wisp_property('barycenter',degree=0).values(n) - topomesh.wisp_property('barycenter',degree=0)[v] for n in vertex_cell_neighbors[v]] for v in topomesh.wisps(0)])
# vertex_cell_laplacian = np.array([[np.sum(vectors,axis=0)/valency for vectors,valency in zip(vertex_cell_edge_vectors[v],vertex_cell_valency[v])] for v in topomesh.wisps(0)])
# vertex_cell_barycenters = np.array([topomesh.wisp_property('barycenter',degree=3).values(vertex_cells[v]) for v in topomesh.wisps(0)])
# vertex_cell_directions = np.array([topomesh.wisp_property('barycenter',degree=0)[v] - vertex_cell_barycenters[v] for v in topomesh.wisps(0)])
# vertex_cell_directions = np.array([vertex_cell_directions[v]/np.linalg.norm(vertex_cell_directions[v],axis=1)[:,np.newaxis] for v in topomesh.wisps(0)])
# vertex_cell_curvature_sign = np.array([-np.sign(np.einsum('ij,ij->i',vertex_cell_laplacian[v],vertex_cell_directions[v])) for v in topomesh.wisps(0)])
# vertex_cell_curvatures = np.array([vertex_cell_curvature_sign[v]*np.linalg.norm(vertex_cell_laplacian[v],axis=1) for v in topomesh.wisps(0)])
# topomesh.update_wisp_property('curvature',degree=degree,values=vertex_cell_curvatures,keys=np.array(list(topomesh.wisps(degree))))
# if property_name == 'gaussian_curvature':
# assert degree==0
# if not 'gaussian_curvature' in topomesh.wisp_property_names(degree):
# topomesh.add_wisp_property('curvature',degree=degree)
# if not topomesh.has_wisp_property('borders',degree=3,is_computed=True):
# compute_topomesh_property(topomesh,'borders',degree=3)
# if not topomesh.has_wisp_property('vertices',degree=3,is_computed=True):
# compute_topomesh_property(topomesh,'vertices',degree=3)
# if not topomesh.has_wisp_property('vertices',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'vertices',degree=2)
# if not topomesh.has_wisp_property('cells',degree=0,is_computed=True):
# compute_topomesh_property(topomesh,'cells',degree=0)
# if not topomesh.has_wisp_property('length',degree=1,is_computed=True):
# compute_topomesh_property(topomesh,'length',degree=1)
# vertex_cells = topomesh.wisp_property('cells',0).values()
# vertex_cell_curvatures = np.array([{} for v in topomesh.wisps(0)])
# cell_triangles = np.array(np.concatenate(topomesh.wisp_property('borders',degree=3).values()),int)
# cell_triangle_cells = np.array(np.concatenate([[c for t in topomesh.wisp_property('borders',degree=3)[c]] for c in topomesh.wisps(3)]),int)
# triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(cell_triangles)
# # rotated_triangle_vertices = np.transpose([triangle_vertices[:,2],triangle_vertices[:,0],triangle_vertices[:,1]])
# # antirotated_triangle_vertices = np.transpose([triangle_vertices[:,1],triangle_vertices[:,2],triangle_vertices[:,0]])
# # triangle_vertices = np.append(np.append(triangle_vertices,rotated_triangle_vertices,axis=0),antirotated_triangle_vertices,axis=0)
# edge_index_list = np.array([[1, 2],[0, 1],[0, 2]])
# triangle_edge_vertices = triangle_vertices[:,edge_index_list]
# triangle_edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,1]) - topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,0])
# #triangle_edge_lengths = np.power(np.sum(np.power(triangle_edge_vectors,2.0),axis=2),0.5)
# triangle_edge_lengths = np.linalg.norm(triangle_edge_vectors,axis=2)
# triangle_edge_directions = triangle_edge_vectors/triangle_edge_lengths[...,np.newaxis]
# triangle_perimeters = np.sum(triangle_edge_lengths,axis=1)
# triangle_areas = np.sqrt((triangle_perimeters/2.0)*(triangle_perimeters/2.0-triangle_edge_lengths[:,0])*(triangle_perimeters/2.0-triangle_edge_lengths[:,1])*(triangle_perimeters/2.0-triangle_edge_lengths[:,2]))
# triangle_cosines = np.zeros_like(triangle_edge_lengths,np.float32)
# triangle_cosines[:,0] = (triangle_edge_lengths[:,1]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2)/(2.0*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2])
# triangle_cosines[:,1] = (triangle_edge_lengths[:,2]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,1]**2)/(2.0*triangle_edge_lengths[:,2]*triangle_edge_lengths[:,0])
# triangle_cosines[:,2] = (triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)/(2.0*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1])
# triangle_angles = np.arccos(triangle_cosines)
# # for v in topomesh.wisps(0):
# # vertex_angle_sum = nd.sum(triangle_angles[np.where(triangle_vertices==v)],cell_triangle_cells[np.where(triangle_vertices==v)[0]],index=topomesh.wisp_property('cells',0)[v])
# # vertex_area = nd.sum(triangle_areas[np.where(triangle_vertices==v)[0]]/3.,cell_triangle_cells[np.where(triangle_vertices==v)[0]],index=topomesh.wisp_property('cells',0)[v])
# # vertex_curvature = (2*np.pi - vertex_angle_sum)/vertex_area
# # vertex_cell_curvatures[v] = vertex_curvature
# #topomesh.wisp_property('curvature',degree=0)[v] = vertex_curvature
# for c in topomesh.wisps(3):
# vertex_angle_sum = nd.sum(triangle_angles[np.where(cell_triangle_cells==c)],triangle_vertices[np.where(cell_triangle_cells==c)],index=topomesh.wisp_property('vertices',degree=3)[c])
# vertex_area = nd.sum(triangle_areas[np.where(cell_triangle_cells==c)[0]][:,np.newaxis]/3.,triangle_vertices[np.where(cell_triangle_cells==c)],index=topomesh.wisp_property('vertices',degree=3)[c])
# vertex_curvature = array_dict(values=(2*np.pi - vertex_angle_sum)/vertex_area,keys=topomesh.wisp_property('vertices',degree=3)[c])
# # print vertex_curvature
# for v in vertex_curvature.keys():
# vertex_cell_curvatures[v][c] = vertex_curvature[v]
# vertex_cell_curvatures = np.array([np.array(vertex_cell_curvatures[v].values()) for v in topomesh.wisps(0)])
# topomesh.update_wisp_property('gaussian_curvature',degree=degree,values=vertex_cell_curvatures,keys=np.array(list(topomesh.wisps(degree))))
# if property_name == 'mean_curvature':
# assert degree==0
# if not 'curvature' in topomesh.wisp_property_names(degree):
# topomesh.add_wisp_property('curvature',degree=degree)
# if not topomesh.has_wisp_property('borders',degree=3,is_computed=True):
# compute_topomesh_property(topomesh,'borders',degree=3)
# if not topomesh.has_wisp_property('vertices',degree=3,is_computed=True):
# compute_topomesh_property(topomesh,'vertices',degree=3)
# if not topomesh.has_wisp_property('vertices',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'vertices',degree=2)
# if not topomesh.has_wisp_property('cells',degree=0,is_computed=True):
# compute_topomesh_property(topomesh,'cells',degree=0)
# if not topomesh.has_wisp_property('length',degree=1,is_computed=True):
# compute_topomesh_property(topomesh,'length',degree=1)
# if not topomesh.has_wisp_property('barycenter',degree=3,is_computed=True):
# compute_topomesh_property(topomesh,'barycenter',degree=3)
# if not topomesh.has_wisp_property('barycenter',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'barycenter',degree=2)
# if not topomesh.has_wisp_property('area',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'area',degree=2)
# if not topomesh.has_wisp_property('normal',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'normal',degree=2)
# triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(list(topomesh.wisps(2)))
# rotated_triangle_vertices = np.transpose([triangle_vertices[:,2],triangle_vertices[:,0],triangle_vertices[:,1]])
# antirotated_triangle_vertices = np.transpose([triangle_vertices[:,1],triangle_vertices[:,2],triangle_vertices[:,0]])
# triangle_vertices = np.append(np.append(triangle_vertices,rotated_triangle_vertices,axis=0),antirotated_triangle_vertices,axis=0)
# edge_index_list = np.array([[1, 2],[0, 1],[0, 2]])
# triangle_edge_vertices = triangle_vertices[:,edge_index_list]
# triangle_edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,1]) - topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,0])
# #triangle_edge_lengths = np.power(np.sum(np.power(triangle_edge_vectors,2.0),axis=2),0.5)
# triangle_edge_lengths = np.linalg.norm(triangle_edge_vectors,axis=2)
# triangle_edge_directions = triangle_edge_vectors/triangle_edge_lengths[...,np.newaxis]
# triangle_perimeters = np.sum(triangle_edge_lengths,axis=1)
# triangle_areas = np.sqrt((triangle_perimeters/2.0)*(triangle_perimeters/2.0-triangle_edge_lengths[:,0])*(triangle_perimeters/2.0-triangle_edge_lengths[:,1])*(triangle_perimeters/2.0-triangle_edge_lengths[:,2]))
# triangle_cosines = np.zeros_like(triangle_edge_lengths,np.float32)
# triangle_cosines[:,0] = (triangle_edge_lengths[:,1]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2)/(2.0*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2])
# triangle_cosines[:,1] = (triangle_edge_lengths[:,2]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,1]**2)/(2.0*triangle_edge_lengths[:,2]*triangle_edge_lengths[:,0])
# triangle_cosines[:,2] = (triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)/(2.0*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1])
# triangle_angles = np.arccos(triangle_cosines)
# triangle_sinuses = np.zeros_like(triangle_edge_lengths,np.float32)
# triangle_sinuses[:,0] = np.sqrt(np.array(1.0 - np.power(triangle_edge_lengths[:,1]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2,2.0)/np.power(2.0*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2],2.0),np.float16))
# triangle_sinuses[:,1] = np.sqrt(np.array(1.0 - np.power(triangle_edge_lengths[:,2]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,1]**2,2.0)/np.power(2.0*triangle_edge_lengths[:,2]*triangle_edge_lengths[:,0],2.0),np.float16))
# triangle_sinuses[:,2] = np.sqrt(np.array(1.0 - np.power(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2,2.0)/np.power(2.0*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1],2.0),np.float16))
# triangle_sinuses[np.where(triangle_sinuses == 0.0)] = 0.001
# triangle_cotangent_vectors = (triangle_cosines/triangle_sinuses)[...,np.newaxis] * triangle_edge_vectors
# # triangle_cotangent_vectors = 1./np.tan(triangle_angles)[...,np.newaxis] * triangle_edge_vectors
# vertex_cotangent_sum = np.transpose([nd.sum(triangle_cotangent_vectors[:,1,0]+triangle_cotangent_vectors[:,2,0],triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0)))),
# nd.sum(triangle_cotangent_vectors[:,1,1]+triangle_cotangent_vectors[:,2,1],triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0)))),
# nd.sum(triangle_cotangent_vectors[:,1,2]+triangle_cotangent_vectors[:,2,2],triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0))))])
# vertex_area = nd.sum(triangle_areas,triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0))))
# # vertex_area = nd.sum(triangle_cotangent_square_lengths[:,1]+triangle_cotangent_square_lengths[:,2], triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0))))/8.
# # vertex_mean_curvature_vectors = vertex_cotangent_sum/(2.*vertex_area[:,np.newaxis])
# vertex_mean_curvature_vectors = (3.*vertex_cotangent_sum)/(4.*vertex_area[:,np.newaxis])
# triangle_cotangent_vectors = 1./np.tan(triangle_angles)[...,np.newaxis] * triangle_edge_vectors
# triangle_cotangent_square_lengths = 1./np.tan(triangle_angles) * triangle_edge_lengths**2
# # triangle_cell_directions = topomesh.wisp_property('barycenter',degree=2).values(list(topomesh.wisps(2))) - topomesh.wisp_property('barycenter',degree=3).values([topomesh.wisp_property('cells',degree=2)[t][0] for t in topomesh.wisps(2)])
# # triangle_cell_directions = np.tile(triangle_cell_directions,(3,1))
# # triangle_normal_vectors = np.cross(triangle_edge_vectors[:,0],triangle_edge_vectors[:,1])
# # triangle_normal_vectors = triangle_normal_vectors/np.linalg.norm(triangle_normal_vectors,axis=1)[:,np.newaxis]
# # reversed_normals = np.where(np.einsum('ij,ij->i',triangle_normal_vectors,triangle_cell_directions) < 0)[0]
# # triangle_normal_vectors[reversed_normals] = -triangle_normal_vectors[reversed_normals]
# triangle_normal_vectors = np.tile(topomesh.wisp_property('normal',degree=2).values(list(topomesh.wisps(2))),(3,1))
# triangle_areas = np.tile(topomesh.wisp_property('area',degree=2).values(list(topomesh.wisps(2))),(3))
# vertex_normal_vectors = np.transpose([nd.sum(triangle_normal_vectors[:,0]*triangle_areas,triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0)))),
# nd.sum(triangle_normal_vectors[:,1]*triangle_areas,triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0)))),
# nd.sum(triangle_normal_vectors[:,2]*triangle_areas,triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0))))])
# vertex_mean_curvature_sign = -np.sign(np.einsum('ij,ij->i',vertex_mean_curvature_vectors,vertex_normal_vectors))
# vertex_mean_curvatures = array_dict(values=vertex_mean_curvature_sign*np.linalg.norm(vertex_mean_curvature_vectors,axis=1),keys=np.array(list(topomesh.wisps(0))))
# topomesh.update_wisp_property('mean_curvature',degree=degree,values=vertex_mean_curvatures,keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'length':
assert degree == 1
if not 'length' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('length',degree=degree)
if not topomesh.has_wisp_property('borders',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'borders',degree=degree)
vertices_positions = positions.values(topomesh.wisp_property('borders',degree).values(list(topomesh.wisps(degree))))
edge_vectors = vertices_positions[:,1] - vertices_positions[:,0]
edge_lengths = np.power(np.sum(np.power(edge_vectors,2.0),axis=1),0.5)
topomesh.update_wisp_property('length',degree=degree,values=edge_lengths,keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'perimeter':
assert degree == 2
if not 'perimeter' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('perimeter',degree=degree)
if not topomesh.has_wisp_property('borders',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'borders',degree=degree)
if not topomesh.has_wisp_property('length',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'length',degree=1)
edge_lengths = topomesh.wisp_property('length',degree=1).values(topomesh.wisp_property('borders',degree=degree).values(list(topomesh.wisps(degree))))
triangle_perimeters = np.sum(edge_lengths,axis=1)
topomesh.update_wisp_property('perimeter',degree=degree,values=triangle_perimeters,keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'area':
assert degree == 2
if not 'area' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('area',degree=degree)
if is_triangular(topomesh):
if not topomesh.has_wisp_property('borders',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'borders',degree=degree)
if not topomesh.has_wisp_property('length',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'length',degree=1)
edge_lengths = topomesh.wisp_property('length',degree=1).values(topomesh.wisp_property('borders',degree=degree).values(list(topomesh.wisps(degree))))
triangle_perimeters = np.sum(edge_lengths,axis=1)
triangle_areas = np.sqrt(np.maximum(0,(triangle_perimeters/2.0)*(triangle_perimeters/2.0-edge_lengths[:,0])*(triangle_perimeters/2.0-edge_lengths[:,1])*(triangle_perimeters/2.0-edge_lengths[:,2])))
topomesh.update_wisp_property('area',degree=degree,values=triangle_areas,keys=np.array(list(topomesh.wisps(degree))))
else:
if not topomesh.has_wisp_property('oriented_vertices',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'oriented_vertices',degree=2)
face_oriented_points = topomesh.wisp_property('barycenter',0).values(topomesh.wisp_property('oriented_vertices',2).values(list(topomesh.wisps(2))))
face_center_points = [np.mean(p,axis=0) for p in face_oriented_points]
face_oriented_triangle_faces = np.concatenate([[f for p in points] for f,points in zip(topomesh.wisps(2),face_oriented_points)])
# face_oriented_triangles = np.concatenate([[[p1,p2,c] for p1,p2 in zip(p,list(p[1:])+[p[0]])] for p,c in zip(face_oriented_points,face_center_points)])
# triangle_edge_list = np.array([[1, 2],[2, 0],[0, 1]])
# face_oriented_triangle_edge_lengths = np.linalg.norm(face_oriented_triangles[:,triangle_edge_list][:,:,1] - face_oriented_triangles[:,triangle_edge_list][:,:,0],axis=2)
# face_oriented_triangle_perimeters = np.sum(face_oriented_triangle_edge_lengths,axis=1)
# face_oriented_triangle_areas = np.sqrt(np.maximum(0, (face_oriented_triangle_perimeters / 2.0) * (face_oriented_triangle_perimeters / 2.0 - face_oriented_triangle_edge_lengths[:, 0]) * (face_oriented_triangle_perimeters / 2.0 - face_oriented_triangle_edge_lengths[:, 1]) * (face_oriented_triangle_perimeters / 2.0 - face_oriented_triangle_edge_lengths[:, 2])))
# face_areas = nd.sum(face_oriented_triangle_areas,face_oriented_triangle_faces,index=list(topomesh.wisps(2)))
face_oriented_center_vectors = np.concatenate([[[p1-c,p2-c] for p1,p2 in zip(p,list(p[1:])+[p[0]])] for p,c in zip(face_oriented_points,face_center_points)])
face_oriented_cross_products = np.linalg.norm(np.cross(face_oriented_center_vectors[:,0],face_oriented_center_vectors[:,1]),axis=1)/2.
face_areas = nd.sum(face_oriented_cross_products,face_oriented_triangle_faces,index=list(topomesh.wisps(2)))
topomesh.update_wisp_property('area', degree=degree, values=face_areas, keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'eccentricity':
assert degree == 2
if not 'eccentricity' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('eccentricity',degree=degree)
if not topomesh.has_wisp_property('vertices',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=2)
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(list(topomesh.wisps(degree)))
edge_index_list = np.array([[1, 2],[0, 1],[0, 2]])
triangle_edge_vertices = triangle_vertices[:,edge_index_list]
triangle_edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,1]) - topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,0])
triangle_edge_lengths = np.linalg.norm(triangle_edge_vectors,axis=2)
triangle_sinuses = np.zeros_like(triangle_edge_lengths,np.float64)
triangle_sinuses[:,0] = np.sqrt(np.maximum(0,np.array(1.0 - np.power((triangle_edge_lengths[:,1]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2)/np.maximum(1e-5,2.0*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2]),2.0),np.float64)))
triangle_sinuses[:,1] = np.sqrt(np.maximum(0,np.array(1.0 - np.power((triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)/np.maximum(1e-5,2.0*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1]),2.0),np.float64)))
triangle_sinuses[:,2] = np.sqrt(np.maximum(0,np.array(1.0 - np.power((triangle_edge_lengths[:,2]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,1]**2)/np.maximum(1e-5,2.0*triangle_edge_lengths[:,2]*triangle_edge_lengths[:,0]),2.0),np.float64)))
triangle_sinus_eccentricities = 1.0 - (2.0*(triangle_sinuses[:,0]+triangle_sinuses[:,1]+triangle_sinuses[:,2]))/(3.*np.sqrt(3.))
topomesh.update_wisp_property('eccentricity',degree=degree,values=triangle_sinus_eccentricities,keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'skewness':
assert degree == 2
if not property_name in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property(property_name=property_name, degree=degree)
if not topomesh.has_wisp_property('length',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'length',degree=1)
if not topomesh.has_wisp_property('borders',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'borders',degree=degree)
edge_lengths = topomesh.wisp_property('length',degree=1).values(topomesh.wisp_property('borders', degree=degree).values(list(topomesh.wisps(degree))))
min_length, max_length = np.min(edge_lengths, axis=1), np.max(edge_lengths, axis=1)
skewness = np.divide(min_length, max_length)
topomesh.update_wisp_property(property_name=property_name, degree=degree, values=skewness, keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'roughness':
assert degree == 0
if not topomesh.has_wisp_property(property_name='area', degree=2):
compute_topomesh_property(topomesh, 'area', degree=2)
if not topomesh.has_wisp_property(property_name, degree):
topomesh.add_wisp_property(property_name=property_name, degree=degree)
areas = topomesh.wisp_property('area', 2)
roughness = {}
for vid in topomesh.wisps(degree):
neighbor_areas = areas.values(list(topomesh.regions(0, vid, 2)))
roughness[vid] = min(neighbor_areas) / max(neighbor_areas)
topomesh.update_wisp_property(property_name=property_name, degree=degree, values=roughness, keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'vertex_neighbors':
assert degree == 2
if not topomesh.has_wisp_property(property_name, degree):
topomesh.add_wisp_property(property_name=property_name, degree=degree)
neighbors = {}
for vid in topomesh.wisps(0):
fids = list(topomesh.regions(0, vid, 2))
for fid in fids:
fids_to_add = np.setdiff1d(fids, fid).tolist()
if fid not in neighbors:
neighbors[fid] = fids_to_add
else:
neighbors[fid] += fids_to_add
vertex_neighbors = {fid: np.unique(fids).tolist() for fid, fids in neighbors.iteritems()}
topomesh.update_wisp_property(property_name=property_name, degree=degree, values=vertex_neighbors, keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'normal':
if degree == 2:
if not 'normal' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('normal',degree=degree)
if not topomesh.has_wisp_property('vertices',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=degree)
if not topomesh.has_wisp_property('epidermis',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',2)
if not topomesh.has_wisp_property('barycenter',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',2)
if not topomesh.has_wisp_property('barycenter',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',3)
if normal_method == "orientation":
normal_cell = max(topomesh.wisps(3))+1
topomesh.add_wisp(3,normal_cell)
for t in topomesh.wisps(2):
topomesh.link(3,normal_cell,t)
compute_topomesh_property(topomesh,'oriented_borders',2)
compute_topomesh_property(topomesh,'oriented_borders',3)
triangle_orientations = array_dict(topomesh.wisp_property('oriented_borders',3)[normal_cell][1], topomesh.wisp_property('oriented_borders',3)[normal_cell][0])
topomesh.update_wisp_property('orientation',2,triangle_orientations.values(list(topomesh.wisps(2))),list(topomesh.wisps(2)))
topomesh.remove_wisp(3,normal_cell)
compute_topomesh_property(topomesh,'vertices',2)
compute_topomesh_property(topomesh,'oriented_vertices',2)
vertices_positions = topomesh.wisp_property('barycenter',0).values(topomesh.wisp_property('oriented_vertices',2).values())
normal_vectors = np.cross(vertices_positions[:,1]-vertices_positions[:,0],vertices_positions[:,2]-vertices_positions[:,0])
normal_norms = np.linalg.norm(normal_vectors,axis=1)
normal_orientations = topomesh.wisp_property('orientation',2).values()
face_vectors = vertices_positions.mean(axis=1) - topomesh.wisp_property('barycenter',0).values().mean(axis=0)
global_normal_orientation = np.sign(np.einsum("...ij,...ij->...i",normal_orientations[:,np.newaxis]*normal_vectors,face_vectors).mean())
global_normal_orientation += (global_normal_orientation==0)
topomesh.update_wisp_property('normal',2,global_normal_orientation*normal_orientations[:,np.newaxis]*normal_vectors/normal_norms[:,np.newaxis],list(topomesh.wisps(2)))
elif normal_method == "density":
vertices_positions = positions.values(topomesh.wisp_property('vertices',degree).values(list(topomesh.wisps(degree))))
normal_vectors = np.cross(vertices_positions[:,1]-vertices_positions[:,0],vertices_positions[:,2]-vertices_positions[:,0])
# reversed_normals = np.where(normal_vectors[:,2] < 0)[0]
from cellcomplex.property_topomesh.utils.implicit_surfaces import point_spherical_density
if object_positions is None:
object_positions = topomesh.wisp_property('barycenter',3)
triangle_epidermis = topomesh.wisp_property('epidermis',2).values()
triangle_exterior_density = point_spherical_density(object_positions,topomesh.wisp_property('barycenter',2).values()[triangle_epidermis]+normal_vectors[triangle_epidermis],sphere_radius=object_radius,k=0.5)
triangle_interior_density = point_spherical_density(object_positions,topomesh.wisp_property('barycenter',2).values()[triangle_epidermis]-normal_vectors[triangle_epidermis],sphere_radius=object_radius,k=0.5)
normal_orientation = 2*(triangle_exterior_density<triangle_interior_density)-1
normal_vectors[triangle_epidermis] = normal_orientation[...,np.newaxis]*normal_vectors[triangle_epidermis]
normal_norms = np.linalg.norm(normal_vectors,axis=1)
# normal_norms[np.where(normal_norms==0)] = 0.001
topomesh.update_wisp_property('normal',degree=degree,values=normal_vectors/np.maximum(1e-5,normal_norms)[:,np.newaxis],keys=np.array(list(topomesh.wisps(degree))))
elif normal_method == "barycenter":
vertices_positions = positions.values(topomesh.wisp_property('vertices',degree).values(list(topomesh.wisps(degree))))
normal_vectors = np.cross(vertices_positions[:,1]-vertices_positions[:,0],vertices_positions[:,2]-vertices_positions[:,0])
barycenter_vectors = topomesh.wisp_property('barycenter',2).values() - positions.values().mean(axis=0)
normal_orientation = np.sign(np.einsum("...ij,...ij->...i",normal_vectors,barycenter_vectors))
normal_vectors = normal_orientation[...,np.newaxis]*normal_vectors
normal_norms = np.linalg.norm(normal_vectors,axis=1)
topomesh.update_wisp_property('normal',degree=degree,values=normal_vectors/normal_norms[:,np.newaxis],keys=np.array(list(topomesh.wisps(degree))))
elif degree == 1:
if not 'normal' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('normal',degree=degree)
if not topomesh.has_wisp_property('vertices',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=1)
if not topomesh.has_wisp_property('regions',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'regions',degree=1)
if not topomesh.has_wisp_property('barycenter',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',1)
if not topomesh.has_wisp_property('barycenter',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',2)
edge_n_faces = np.array([topomesh.nb_regions(1,e) for e in topomesh.wisps(1)])
edge_points = topomesh.wisp_property('barycenter',0).values(topomesh.wisp_property('vertices',1).values(list(topomesh.wisps(1))))
edge_vectors = edge_points[:,1] - edge_points[:,0]
edge_face_centers = topomesh.wisp_property('barycenter',2).values(topomesh.wisp_property('regions',1).values(list(topomesh.wisps(1))))
edge_face_centers = np.array(list(map(lambda p : np.mean(p,axis=0), edge_face_centers)))
edge_centers = topomesh.wisp_property('barycenter',1).values(list(topomesh.wisps(1)))
edge_center_vectors = edge_face_centers - edge_centers
edge_plane_normals = np.cross(edge_vectors,edge_center_vectors)
edge_normals = np.cross(edge_vectors,edge_plane_normals)
edge_normals = edge_normals/np.linalg.norm(edge_normals,axis=1)[:,np.newaxis]
edge_normal_orientation = -np.sign(np.einsum("...ij,...ij->...i",edge_normals,edge_center_vectors))
edge_normal_orientation[edge_normal_orientation==0] = -1
edge_normals = edge_normal_orientation[:,np.newaxis] * (edge_n_faces==1)[:,np.newaxis] * edge_normals
topomesh.update_wisp_property('normal', degree=degree, values=edge_normals, keys=np.array(list(topomesh.wisps(degree))))
elif degree == 0:
if not topomesh.has_wisp_property('normal',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'normal',degree=2)
if not topomesh.has_wisp_property('area',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'area',degree=2)
if not topomesh.has_wisp_property('vertices',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=2)
if not topomesh.has_wisp_property('epidermis',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',degree=2)
epidermis_triangles = np.array(list(topomesh.wisps(2)))[topomesh.wisp_property('epidermis',2).values(list(topomesh.wisps(2)))]
# print epidermis_triangles
#triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(list(topomesh.wisps(2)))
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(epidermis_triangles)
rotated_triangle_vertices = np.transpose([triangle_vertices[:,2],triangle_vertices[:,0],triangle_vertices[:,1]])
antirotated_triangle_vertices = np.transpose([triangle_vertices[:,1],triangle_vertices[:,2],triangle_vertices[:,0]])
triangle_vertices = np.append(np.append(triangle_vertices,rotated_triangle_vertices,axis=0),antirotated_triangle_vertices,axis=0)
# triangle_areas = np.tile(topomesh.wisp_property('area',degree=2).values(list(topomesh.wisps(2))),(3))
# triangle_normal_vectors = np.tile(topomesh.wisp_property('normal',degree=2).values(list(topomesh.wisps(2))),(3,1))
triangle_areas = np.tile(topomesh.wisp_property('area',degree=2).values(epidermis_triangles),(3))
triangle_normal_vectors = np.tile(topomesh.wisp_property('normal',degree=2).values(epidermis_triangles),(3,1))
vertex_normal_vectors = np.transpose([nd.sum(triangle_normal_vectors[:,k]*triangle_areas,triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0)))) for k in [0,1,2]])
normal_norms = np.linalg.norm(vertex_normal_vectors,axis=1)
normal_norms[np.where(normal_norms==0)] = 0.001
topomesh.update_wisp_property('normal',degree=degree,values=vertex_normal_vectors/normal_norms[:,np.newaxis],keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'angles':
assert degree == 2
if not 'angles' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('angles',degree=degree)
if not topomesh.has_wisp_property('vertices',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=2)
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(list(topomesh.wisps(degree)))
edge_index_list = np.array([[1, 2],[0, 1],[0, 2]])
triangle_edge_vertices = triangle_vertices[:,edge_index_list]
triangle_edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,1]) - topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,0])
triangle_edge_lengths = np.linalg.norm(triangle_edge_vectors,axis=2)
triangle_cosines = np.zeros_like(triangle_edge_lengths,np.float32)
triangle_cosines[:,0] = (triangle_edge_lengths[:,1]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2)/np.maximum(1e-5,2.0*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2])
triangle_cosines[:,1] = (triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)/np.maximum(1e-5,2.0*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1])
triangle_cosines[:,2] = (triangle_edge_lengths[:,2]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,1]**2)/np.maximum(1e-5,2.0*triangle_edge_lengths[:,2]*triangle_edge_lengths[:,0])
triangle_angles = np.arccos(triangle_cosines)
topomesh.update_wisp_property('angles',degree=degree,values=triangle_angles,keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'incircle center':
assert degree == 2
if not topomesh.has_wisp_property('vertices',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=2)
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(list(topomesh.wisps(degree)))
edge_index_list = np.array([[1, 2],[0, 2],[0, 1]])
triangle_edge_vertices = triangle_vertices[:,edge_index_list]
triangle_edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,1]) - topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,0])
triangle_edge_lengths = np.linalg.norm(triangle_edge_vectors,axis=2)
triangle_sinuses = np.zeros_like(triangle_edge_lengths,np.float32)
triangle_sinuses[:,0] = np.sqrt(1.0 - np.power((triangle_edge_lengths[:,1]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2)/(2.0*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2]),2.0))
triangle_sinuses[:,1] = np.sqrt(1.0 - np.power((triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)/(2.0*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1]),2.0))
triangle_sinuses[:,2] = np.sqrt(1.0 - np.power((triangle_edge_lengths[:,2]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,1]**2)/(2.0*triangle_edge_lengths[:,2]*triangle_edge_lengths[:,0]),2.0))
# print triangle_sinuses.shape
# print triangle_vertices.shape
triangle_centers = ((topomesh.wisp_property('barycenter',degree=0).values(triangle_vertices)/triangle_sinuses[...,np.newaxis]).sum(axis=1))/((1./triangle_sinuses).sum(axis=1)[...,np.newaxis])
# print triangle_centers
# print triangle_centers.shape
topomesh.update_wisp_property('incircle center',degree=degree,values=triangle_centers,keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'volume':
assert degree == 3
if not 'volume' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('volume',degree=degree)
if not topomesh.has_wisp_property('borders',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'borders',degree=degree)
if not topomesh.has_wisp_property('vertices',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=2)
if not topomesh.has_wisp_property('barycenter',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=degree)
cell_triangles = np.array(np.concatenate(topomesh.wisp_property('borders',degree=3).values(list(topomesh.wisps(degree)))),int)
cell_tetrahedra_cells = np.array(np.concatenate([[c for t in topomesh.wisp_property('borders',degree=3)[c]] for c in topomesh.wisps(3)]),int)
cell_triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(cell_triangles)
cell_tetrahedra = np.concatenate([topomesh.wisp_property('barycenter',degree=3).values(cell_tetrahedra_cells)[:,np.newaxis], topomesh.wisp_property('barycenter',degree=0).values(cell_triangle_vertices)],axis=1)
cell_tetra_matrix = np.transpose(np.array([cell_tetrahedra[:,1],cell_tetrahedra[:,2],cell_tetrahedra[:,3]]) - cell_tetrahedra[:,0],axes=(1,2,0))
cell_tetra_volume = abs(np.linalg.det(cell_tetra_matrix))/6.0
cell_volumes = nd.sum(cell_tetra_volume,cell_tetrahedra_cells,index=list(topomesh.wisps(3)))
topomesh.update_wisp_property('volume',degree=degree,values=cell_volumes,keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'convexhull_volume':
import openalea.plantgl.all as pgl
assert degree == 3
if not 'convexhull_volume' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('convexhull_volume',degree=degree)
if not topomesh.has_wisp_property('vertices',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=degree)
for c in topomesh.wisps(3):
cell_vertices = topomesh.wisp_property('vertices',degree)[c]
if len(cell_vertices)>0:
convex_hull = pgl.Fit(pgl.Point3Array(positions.values(cell_vertices))).convexHull()
topomesh.wisp_property('convexhull_volume',degree=degree)[c] = pgl.volume(convex_hull)
if property_name == 'surface':
assert degree == 3
if not 'surface' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('surface',degree=degree)
if not topomesh.has_wisp_property('borders',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'borders',degree=degree)
if not topomesh.has_wisp_property('area',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'area',degree=2)
cell_triangles = np.array(np.concatenate(topomesh.wisp_property('borders',degree=3).values(list(topomesh.wisps(degree)))),int)
cell_triangle_cells = np.array(np.concatenate([[c for t in topomesh.wisp_property('borders',degree=3)[c]] for c in topomesh.wisps(3)]),int)
cell_triangle_areas = topomesh.wisp_property('area',degree=2).values(cell_triangles)
cell_surfaces = nd.sum(cell_triangle_areas,cell_triangle_cells,index=list(topomesh.wisps(3)))
topomesh.update_wisp_property('surface',degree=degree,values=cell_surfaces,keys=np.array(list(topomesh.wisps(degree))))
if property_name == 'convexhull_surface':
import openalea.plantgl.all as pgl
assert degree == 3
if not 'convexhull_surface' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('convexhull_surface',degree=degree)
if not topomesh.has_wisp_property('vertices',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=degree)
for c in topomesh.wisps(3):
cell_vertices = topomesh.wisp_property('vertices',degree)[c]
if len(cell_vertices)>0:
convex_hull = pgl.Fit(pgl.Point3Array(topomesh.wisp_property('barycenter',0).values(cell_vertices))).convexHull()
topomesh.wisp_property('convexhull_surface',degree=degree)[c] = pgl.surface(convex_hull)
if property_name == 'convexity':
assert degree == 3
if not 'convexity' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('convexity',degree=degree)
if not topomesh.has_wisp_property('volume',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'volume',degree=degree)
if not topomesh.has_wisp_property('convexhull_volume',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'convexhull_volume',degree=degree)
topomesh.update_wisp_property('convexity',degree=degree,values=topomesh.wisp_property('volume',3).values(topomesh.wisp_property('convexhull_volume',3).keys())/topomesh.wisp_property('convexhull_volume',3).values(topomesh.wisp_property('convexhull_volume',3).keys()),keys=topomesh.wisp_property('convexhull_volume',3).keys())
if property_name == 'epidermis':
if not 'epidermis' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('epidermis',degree=degree)
if degree != 2:
if not topomesh.has_wisp_property('faces',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'faces',degree=degree)
if not topomesh.has_wisp_property('epidermis',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',degree=2)
epidermis = np.array([topomesh.wisp_property('epidermis',degree=2).values(topomesh.wisp_property('faces',degree=degree)[w]).any() for w in topomesh.wisps(degree)])
topomesh.update_wisp_property('epidermis',degree=degree,values=epidermis,keys=np.array(list(topomesh.wisps(degree))))
else:
if not topomesh.has_wisp_property('cells',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'cells',degree=2)
topomesh.update_wisp_property('epidermis',degree=2,values=(np.array(list(map(len,topomesh.wisp_property('cells',2).values(list(topomesh.wisps(2)))))) == 1),keys=np.array(list(topomesh.wisps(2))))
# if property_name == 'epidermis_curvature':
# if degree == 0:
# if not topomesh.has_wisp_property('epidermis',degree=0,is_computed=True):
# compute_topomesh_property(topomesh,'epidermis',degree=0)
# if not topomesh.has_wisp_property('vertices',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'vertices',degree=2)
# if not topomesh.has_wisp_property('barycenter',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'barycenter',degree=2)
# if not topomesh.has_wisp_property('barycenter',degree=3,is_computed=True):
# compute_topomesh_property(topomesh,'barycenter',degree=3)
# if not topomesh.has_wisp_property('area',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'area',degree=2)
# if not topomesh.has_wisp_property('normal',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'normal',degree=2)
# curvature_radius = kwargs.get('radius',15)
# epidermis_vertices = np.array(list(topomesh.wisps(0)))[topomesh.wisp_property('epidermis',0).values(list(topomesh.wisps(0)))]
# vertex_distances = np.array([vq(topomesh.wisp_property('barycenter',degree=0).values(epidermis_vertices),topomesh.wisp_property('barycenter',degree=0).values([v]))[1] for v in epidermis_vertices])
# vertex_neighborhood = np.array([np.array(list(topomesh.wisps(0)))[np.where(vertex_distances[v]<curvature_radius)] for v in range(epidermis_vertices.shape[0])])
# vertex_neighborhood_barycenter = np.array([np.mean(topomesh.wisp_property('barycenter',degree=0).values(vertex_neighborhood[v]),axis=0) for v in range(epidermis_vertices.shape[0])])
# vertex_vectors = topomesh.wisp_property('barycenter',degree=0).values(epidermis_vertices) - vertex_neighborhood_barycenter
# vertex_vector_norm = np.linalg.norm(vertex_vectors,axis=1)
# epidermis_triangles = np.array(list(topomesh.wisps(2)))[topomesh.wisp_property('epidermis',2).values(list(topomesh.wisps(2)))]
# triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(epidermis_triangles)
# rotated_triangle_vertices = np.transpose([triangle_vertices[:,2],triangle_vertices[:,0],triangle_vertices[:,1]])
# antirotated_triangle_vertices = np.transpose([triangle_vertices[:,1],triangle_vertices[:,2],triangle_vertices[:,0]])
# triangle_vertices = np.append(np.append(triangle_vertices,rotated_triangle_vertices,axis=0),antirotated_triangle_vertices,axis=0)
# triangle_normal_vectors = np.tile(topomesh.wisp_property('normal',degree=2).values(epidermis_triangles),(3,1))
# triangle_areas = np.tile(topomesh.wisp_property('area',degree=2).values(epidermis_triangles),(3))
# print triangle_normal_vectors.shape
# print triangle_areas.shape
# print triangle_vertices.shape
# vertex_normal_vectors = np.transpose([nd.sum(triangle_normal_vectors[:,0]*triangle_areas,triangle_vertices[:,0],index=epidermis_vertices),
# nd.sum(triangle_normal_vectors[:,1]*triangle_areas,triangle_vertices[:,0],index=epidermis_vertices),
# nd.sum(triangle_normal_vectors[:,2]*triangle_areas,triangle_vertices[:,0],index=epidermis_vertices)])
# vertex_curvature_sign = np.sign(np.einsum('ij,ij->i',vertex_vectors,vertex_normal_vectors))
# vertex_curvature = array_dict(4.*vertex_curvature_sign*vertex_vector_norm/curvature_radius,keys=epidermis_vertices)
# vertex_curvature_inner = np.array([vertex_curvature[v] if v in epidermis_vertices else 0.0 for v in topomesh.wisps(0)])
# topomesh.update_wisp_property('epidermis_curvature',0,values=vertex_curvature_inner,keys=np.array(list(topomesh.wisps(0))))
# if degree == 2:
# compute_topomesh_property(topomesh,'epidermis_curvature',degree=0,radius=kwargs.get('radius',60))
# triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(list(topomesh.wisps(2)))
# triangle_curvature = topomesh.wisp_property('epidermis_curvature',degree=0).values(triangle_vertices).mean(axis=1)
# topomesh.update_wisp_property('epidermis_curvature',2,values=triangle_curvature,keys=np.array(list(topomesh.wisps(2))))
if property_name in ['principal_curvatures','principal_curvature_tensor','principal_direction_min','principal_direction_max','principal_curvature_min','principal_curvature_max','mean_curvature','gaussian_curvature']:
if degree == 0:
if not topomesh.has_wisp_property('barycenter',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=2)
if not topomesh.has_wisp_property('barycenter',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=3)
if not topomesh.has_wisp_property('area',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'area',degree=2)
if not topomesh.has_wisp_property('normal',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'normal',degree=2)
if not topomesh.has_wisp_property('normal',degree=0,is_computed=True):
compute_topomesh_property(topomesh,'normal',degree=0)
if not topomesh.has_wisp_property('neighors',degree=0,is_computed=True):
compute_topomesh_property(topomesh,'neighbors',0)
if not topomesh.has_wisp_property('epidermis',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',2)
if not topomesh.has_wisp_property('epidermis',degree=0,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',0)
vertex_normal_vectors = topomesh.wisp_property('normal',0)
epidermis_vertices = np.array(list(topomesh.wisps(0)))[topomesh.wisp_property('epidermis',0).values(list(topomesh.wisps(0)))]
vertex_neighbor_vertices = np.concatenate([[n for n in topomesh.wisp_property('neighbors',0)[v] if n in epidermis_vertices] for v in epidermis_vertices])
vertex_neighbor_vertex = np.concatenate([[v for n in topomesh.wisp_property('neighbors',0)[v] if n in epidermis_vertices] for v in epidermis_vertices])
vertex_neighbor_vectors = topomesh.wisp_property('barycenter',0).values(vertex_neighbor_vertex) - topomesh.wisp_property('barycenter',0).values(vertex_neighbor_vertices)
vertex_neighbor_projected_vectors = vertex_neighbor_vectors - np.einsum('ij,ij->i',vertex_neighbor_vectors,vertex_normal_vectors.values(vertex_neighbor_vertex))[:,np.newaxis]*vertex_normal_vectors.values(vertex_neighbor_vertex)
vertex_neighbor_projected_vectors = vertex_neighbor_projected_vectors/np.linalg.norm(vertex_neighbor_projected_vectors,axis=1)[:,np.newaxis]
vertex_neighbor_directional_curvature = 2.*np.einsum('ij,ij->i',vertex_normal_vectors.values(vertex_neighbor_vertex),vertex_neighbor_vectors)/np.power(np.linalg.norm(vertex_neighbor_vectors,axis=1),2.0)
# epidermis_triangles = np.array(list(topomesh.wisps(2)))[topomesh.wisp_property('epidermis',2).values(list(topomesh.wisps(2)))]
compute_topomesh_property(topomesh,'faces',0)
vertex_neighbor_triangles = np.array([[t for t in topomesh.wisp_property('faces',0)[v] if (t in topomesh.wisp_property('faces',0)[n]) and (topomesh.wisp_property('epidermis',2)[t])] for (v,n) in zip(vertex_neighbor_vertex,vertex_neighbor_vertices)])
vertex_neighbor_triangle_area = topomesh.wisp_property('area',2).values(vertex_neighbor_triangles)
vertex_neighbor_triangle_area = np.array([a.sum() for a in vertex_neighbor_triangle_area])
vertex_neighbor_projected_matrix = np.einsum('...ji,...ij->...ij',vertex_neighbor_projected_vectors[...,np.newaxis],vertex_neighbor_projected_vectors[...,np.newaxis])
vertex_neighbor_curvature_matrix = vertex_neighbor_triangle_area[:,np.newaxis,np.newaxis]*vertex_neighbor_directional_curvature[:,np.newaxis,np.newaxis]*vertex_neighbor_projected_matrix
vertex_curvature_matrix = np.zeros((epidermis_vertices.shape[0],3,3))
for i in range(3):
for j in range(3):
vertex_curvature_matrix[:,i,j] = nd.sum(vertex_neighbor_curvature_matrix[:,i,j],vertex_neighbor_vertex,index=epidermis_vertices)
vertex_curvature_matrix = vertex_curvature_matrix / (nd.sum(vertex_neighbor_triangle_area,vertex_neighbor_vertex,index=epidermis_vertices)[:,np.newaxis,np.newaxis]+np.power(10.,-10))
vertex_curvature_tensor = array_dict(vertex_curvature_matrix,epidermis_vertices)
try:
vertex_curvature_matrix_eigenvalues, vertex_curvature_matrix_eigenvectors = np.linalg.eig(vertex_curvature_matrix)
except np.linalg.LinAlgError:
#except LinAlgError:
ok_curvature_matrix = np.unique(np.where(1-np.isnan(vertex_curvature_matrix))[0])
nan_curvature_matrix = np.unique(np.where(np.isnan(vertex_curvature_matrix))[0])
vertex_curvature_matrix_trunc = np.delete(vertex_curvature_matrix,nan_curvature_matrix)
vertex_curvature_matrix_eigenvalues_trunc, vertex_curvature_matrix_eigenvectors_trunc = np.linalg.eig(vertex_curvature_matrix_trunc)
vertex_curvature_matrix_eigenvalues = np.zeros((vertex_curvature_matrix.shape[0],3))
vertex_curvature_matrix_eigenvectors = np.zeros((vertex_curvature_matrix.shape[0],3,3))
vertex_curvature_matrix_eigenvalues[ok_curvature_matrix] = vertex_curvature_matrix_eigenvalues_trunc
vertex_curvature_matrix_eigenvectors[ok_curvature_matrix] = vertex_curvature_matrix_eigenvectors_trunc
# vertex_principal_curvature_min = array_dict(vertex_curvature_matrix_eigenvalues[tuple([np.arange(len(epidermis_vertices)),np.argsort(np.abs(vertex_curvature_matrix_eigenvalues))[:,1]])],epidermis_vertices)
# vertex_principal_curvature_max = array_dict(vertex_curvature_matrix_eigenvalues[tuple([np.arange(len(epidermis_vertices)),np.argsort(np.abs(vertex_curvature_matrix_eigenvalues))[:,2]])],epidermis_vertices)
# vertex_principal_direction_min = array_dict(vertex_curvature_matrix_eigenvectors[tuple([np.arange(len(epidermis_vertices)),np.argsort(np.abs(vertex_curvature_matrix_eigenvalues))[:,1]])],epidermis_vertices)
# vertex_principal_direction_max = array_dict(vertex_curvature_matrix_eigenvectors[tuple([np.arange(len(epidermis_vertices)),np.argsort(np.abs(vertex_curvature_matrix_eigenvalues))[:,2]])],epidermis_vertices)
vertex_principal_curvature_min = array_dict([np.sort(2*v[np.argsort(np.abs(v))[1:]])[0] for v in vertex_curvature_matrix_eigenvalues],epidermis_vertices)
vertex_principal_curvature_max = array_dict([np.sort(2*v[np.argsort(np.abs(v))[1:]])[1] for v in vertex_curvature_matrix_eigenvalues],epidermis_vertices)
vertex_principal_direction_min = array_dict([vec[np.argsort(v[np.argsort(np.abs(v))[1:]])[0]] for vec,v in zip(vertex_curvature_matrix_eigenvectors,vertex_curvature_matrix_eigenvalues)],epidermis_vertices)
vertex_principal_direction_max = array_dict([vec[np.argsort(v[np.argsort(np.abs(v))[1:]])[1]] for vec,v in zip(vertex_curvature_matrix_eigenvectors,vertex_curvature_matrix_eigenvalues)],epidermis_vertices)
topomesh.update_wisp_property('principal_curvature_tensor',0,np.array([vertex_curvature_tensor[v] if v in epidermis_vertices else np.zeros((3,3),float) for v in topomesh.wisps(0)]),np.array(list(topomesh.wisps(0))))
topomesh.update_wisp_property('principal_direction_min',0,np.array([vertex_principal_direction_min[v] if v in epidermis_vertices else np.array([0.,0.,0.]) for v in topomesh.wisps(0)]),np.array(list(topomesh.wisps(0))))
topomesh.update_wisp_property('principal_direction_max',0,np.array([vertex_principal_direction_max[v] if v in epidermis_vertices else np.array([0.,0.,0.]) for v in topomesh.wisps(0)]),np.array(list(topomesh.wisps(0))))
topomesh.update_wisp_property('principal_curvature_min',0,np.array([vertex_principal_curvature_min[v] if v in epidermis_vertices else 0. for v in topomesh.wisps(0)]),np.array(list(topomesh.wisps(0))))
topomesh.update_wisp_property('principal_curvature_max',0,np.array([vertex_principal_curvature_max[v] if v in epidermis_vertices else 0. for v in topomesh.wisps(0)]),np.array(list(topomesh.wisps(0))))
topomesh.update_wisp_property('mean_curvature',0,(topomesh.wisp_property('principal_curvature_max',0).values(list(topomesh.wisps(0)))+topomesh.wisp_property('principal_curvature_min',0).values(list(topomesh.wisps(0))))/2.,np.array(list(topomesh.wisps(0))))
topomesh.update_wisp_property('gaussian_curvature',0,topomesh.wisp_property('principal_curvature_max',0).values(list(topomesh.wisps(0)))*topomesh.wisp_property('principal_curvature_min',0).values(list(topomesh.wisps(0))),np.array(list(topomesh.wisps(0))))
elif degree == 2:
# if not topomesh.has_wisp_property('principal_curvature_tensor',degree=0,is_computed=True):
# compute_topomesh_property(topomesh,'principal_curvature_tensor',degree=0)
# if not topomesh.has_wisp_property('epidermis',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'epidermis',2)
# if not topomesh.has_wisp_property('vertices',degree=2,is_computed=True):
# compute_topomesh_property(topomesh,'vertices',2)
# epidermis_triangles = np.array(list(topomesh.wisps(2)))[topomesh.wisp_property('epidermis',2).values(list(topomesh.wisps(2)))]
# triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(epidermis_triangles)
# triangle_curvature_tensor = array_dict(topomesh.wisp_property('principal_curvature_tensor',degree=0).values(triangle_vertices).mean(axis=1),epidermis_triangles)
# topomesh.update_wisp_property('principal_curvature_tensor',2,np.array([triangle_curvature_tensor[t] if t in epidermis_triangles else np.zeros((3,3),float) for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
# for direction in ['principal_direction_min','principal_direction_max']:
# triangle_direction = array_dict(topomesh.wisp_property(direction,degree=0).values(triangle_vertices).mean(axis=1),epidermis_triangles)
# topomesh.update_wisp_property(direction,2,np.array([triangle_direction[t] if t in epidermis_triangles else np.array([0.,0.,0.]) for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
# for curvature in ['principal_curvature_min','principal_curvature_max','mean_curvature','gaussian_curvature']:
# triangle_curvature = array_dict(topomesh.wisp_property(curvature,degree=0).values(triangle_vertices).mean(axis=1),epidermis_triangles)
# topomesh.update_wisp_property(curvature,2,np.array([triangle_curvature[t] if t in epidermis_triangles else 0. for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
if not topomesh.has_wisp_property('barycenter',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=2)
if not topomesh.has_wisp_property('normal',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'normal',degree=2)
if not topomesh.has_wisp_property('area',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'area',degree=2)
if not topomesh.has_wisp_property('normal',degree=0,is_computed=True):
compute_topomesh_property(topomesh,'normal',degree=0)
epidermis_triangles = np.array(list(topomesh.wisps(2)))[topomesh.wisp_property('epidermis',2).values(list(topomesh.wisps(2))).astype(bool)]
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(epidermis_triangles)
triangle_edge_list = np.array([[1, 2],[0, 2],[0, 1]])
triangle_edge_vertices = triangle_vertices[:,triangle_edge_list]
triangle_edge_points = positions.values(triangle_edge_vertices)
triangle_edge_vectors = triangle_edge_points[:,:,1]-triangle_edge_points[:,:,0]
triangle_vertex_normals = topomesh.wisp_property('normal',0).values(triangle_vertices)
triangle_barycenter_normals = triangle_vertex_normals.mean(axis=1)
#triangle_barycenter_normals = triangle_barycenter_normals/np.linalg.norm(triangle_barycenter_normals,axis=1)[:,np.newaxis]
triangle_barycenter_normal_derivatives = triangle_vertex_normals[:,triangle_edge_list]
triangle_barycenter_normal_derivatives = triangle_barycenter_normal_derivatives[:,:,1] - triangle_barycenter_normal_derivatives[:,:,0]
triangle_barycenter_normal_derivatives = triangle_barycenter_normal_derivatives/np.linalg.norm(triangle_barycenter_normals,axis=1)[:,np.newaxis,np.newaxis]
triangle_barycenter_derivatives_projectors = np.transpose([np.einsum("...ij,...ij->...i",triangle_barycenter_normals,triangle_edge_vectors[:,k])[:,np.newaxis]*triangle_barycenter_normals for k in range(3)],(1,0,2))
triangle_projected_barycenter_derivatives = triangle_edge_vectors - triangle_barycenter_derivatives_projectors
E = np.einsum("...ij,...ij->...i",triangle_projected_barycenter_derivatives[:,1],triangle_projected_barycenter_derivatives[:,1])
F = np.einsum("...ij,...ij->...i",triangle_projected_barycenter_derivatives[:,1],triangle_projected_barycenter_derivatives[:,2])
G = np.einsum("...ij,...ij->...i",triangle_projected_barycenter_derivatives[:,2],triangle_projected_barycenter_derivatives[:,2])
L = -np.einsum("...ij,...ij->...i",triangle_barycenter_normal_derivatives[:,1],triangle_projected_barycenter_derivatives[:,1])
M1 = -np.einsum("...ij,...ij->...i",triangle_barycenter_normal_derivatives[:,1],triangle_projected_barycenter_derivatives[:,2])
M2 = -np.einsum("...ij,...ij->...i",triangle_barycenter_normal_derivatives[:,2],triangle_projected_barycenter_derivatives[:,1])
N = -np.einsum("...ij,...ij->...i",triangle_barycenter_normal_derivatives[:,2],triangle_projected_barycenter_derivatives[:,2])
weingarten_curvature_matrix = np.zeros((len(epidermis_triangles),2,2))
weingarten_curvature_matrix[:,0,0] = (L*G-M1*F)/(E*G-F*F)
weingarten_curvature_matrix[:,0,1] = (M2*G-N*F)/(E*G-F*F)
# weingarten_curvature_matrix[:,0,1] = (M1*E-L*F)/(E*G-F*F)
weingarten_curvature_matrix[:,1,0] = (M1*E-L*F)/(E*G-F*F)
# weingarten_curvature_matrix[:,1,0] = (M2*G-N*F)/(E*G-F*F)
weingarten_curvature_matrix[:,1,1] = (N*E-M2*F)/(E*G-F*F)
weingarten_curvature_matrix_eigenvalues, weingarten_curvature_matrix_eigenvectors = np.linalg.eig(weingarten_curvature_matrix)
weingarten_curvature_matrix_eigenvalues = -weingarten_curvature_matrix_eigenvalues
weingarten_curvature_matrix_eigenvectors = np.transpose(weingarten_curvature_matrix_eigenvectors,(0,2,1))
weingarten_principal_curvature_min = weingarten_curvature_matrix_eigenvalues[tuple([np.arange(len(epidermis_triangles)),np.argsort(np.abs(weingarten_curvature_matrix_eigenvalues))[:,0]])].astype(float)
weingarten_principal_curvature_max = weingarten_curvature_matrix_eigenvalues[tuple([np.arange(len(epidermis_triangles)),np.argsort(np.abs(weingarten_curvature_matrix_eigenvalues))[:,1]])].astype(float)
weingarten_principal_vector_min = weingarten_curvature_matrix_eigenvectors[tuple([np.arange(len(epidermis_triangles)),np.argsort(np.abs(weingarten_curvature_matrix_eigenvalues))[:,0]])].astype(float)
weingarten_principal_vector_max = weingarten_curvature_matrix_eigenvectors[tuple([np.arange(len(epidermis_triangles)),np.argsort(np.abs(weingarten_curvature_matrix_eigenvalues))[:,1]])].astype(float)
# weingarten_principal_curvature_min = weingarten_curvature_matrix_eigenvalues[tuple([np.arange(len(epidermis_triangles)),np.argsort(weingarten_curvature_matrix_eigenvalues)[:,0]])].astype(float)
# weingarten_principal_curvature_max = weingarten_curvature_matrix_eigenvalues[tuple([np.arange(len(epidermis_triangles)),np.argsort(weingarten_curvature_matrix_eigenvalues)[:,1]])].astype(float)
# weingarten_principal_vector_min = weingarten_curvature_matrix_eigenvectors[tuple([np.arange(len(epidermis_triangles)),np.argsort(weingarten_curvature_matrix_eigenvalues)[:,0]])].astype(float)
# weingarten_principal_vector_max = weingarten_curvature_matrix_eigenvectors[tuple([np.arange(len(epidermis_triangles)),np.argsort(weingarten_curvature_matrix_eigenvalues)[:,1]])].astype(float)
weingarten_principal_direction_min = array_dict((weingarten_principal_vector_min[:,:,np.newaxis]*triangle_projected_barycenter_derivatives[:,1:3]).sum(axis=1),epidermis_triangles)
weingarten_principal_direction_max = array_dict((weingarten_principal_vector_max[:,:,np.newaxis]*triangle_projected_barycenter_derivatives[:,1:3]).sum(axis=1),epidermis_triangles)
P = np.transpose([weingarten_principal_direction_max.values(), weingarten_principal_direction_min.values(), triangle_barycenter_normals],(1,2,0))
D = np.array([np.diag(d) for d in np.transpose([weingarten_principal_curvature_max, weingarten_principal_curvature_min, np.zeros_like(epidermis_triangles)])])
P_i = np.array([np.linalg.pinv(p) for p in P])
face_curvature_tensor = np.einsum('...ij,...jk->...ik',P,np.einsum('...ij,...jk->...ik',D,P_i))
face_curvature_matrix_eigenvalues, face_curvature_matrix_eigenvectors = np.linalg.eig(face_curvature_tensor)
# face_principal_curvature_min = array_dict(face_curvature_matrix_eigenvalues[tuple([np.arange(len(epidermis_triangles)),np.argsort(np.abs(face_curvature_matrix_eigenvalues))[:,1]])].astype(float),epidermis_triangles)
# face_principal_curvature_max = array_dict(face_curvature_matrix_eigenvalues[tuple([np.arange(len(epidermis_triangles)),np.argsort(np.abs(face_curvature_matrix_eigenvalues))[:,2]])].astype(float),epidermis_triangles)
face_principal_curvature_min = array_dict(np.array([np.sort(v[np.argsort(np.abs(v))[1:]])[0] for v in face_curvature_matrix_eigenvalues]).astype(float),epidermis_triangles)
face_principal_curvature_max = array_dict(np.array([np.sort(v[np.argsort(np.abs(v))[1:]])[1] for v in face_curvature_matrix_eigenvalues]).astype(float),epidermis_triangles)
# face_curvature_matrix_eigenvectors = np.transpose(face_curvature_matrix_eigenvectors,(0,2,1))
face_curvature_matrix_eigenvectors = np.transpose(face_curvature_matrix_eigenvectors,(0,1,2))
# face_principal_direction_min = array_dict(face_curvature_matrix_eigenvectors[tuple([np.arange(len(epidermis_triangles)),np.argsort(np.abs(face_curvature_matrix_eigenvalues))[:,1]])].astype(float),epidermis_triangles)
# face_principal_direction_max = array_dict(face_curvature_matrix_eigenvectors[tuple([np.arange(len(epidermis_triangles)),np.argsort(np.abs(face_curvature_matrix_eigenvalues))[:,2]])].astype(float),epidermis_triangles)
face_principal_direction_min = array_dict(np.array([vec[np.argsort(v[np.argsort(np.abs(v))[1:]])[0]] for vec,v in zip(face_curvature_matrix_eigenvectors,face_curvature_matrix_eigenvalues)]).astype(float),epidermis_triangles)
face_principal_direction_max = array_dict(np.array([vec[np.argsort(v[np.argsort(np.abs(v))[1:]])[1]] for vec,v in zip(face_curvature_matrix_eigenvectors,face_curvature_matrix_eigenvalues)]).astype(float),epidermis_triangles)
face_principal_curvature_tensor = array_dict(face_curvature_tensor,epidermis_triangles)
topomesh.update_wisp_property('principal_curvature_tensor',2,np.array([face_principal_curvature_tensor[t] if t in epidermis_triangles else np.zeros((3,3)) for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
topomesh.update_wisp_property('principal_direction_min',2,np.array([face_principal_direction_min[t] if t in epidermis_triangles else np.zeros(3) for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
topomesh.update_wisp_property('principal_direction_max',2,np.array([face_principal_direction_max[t] if t in epidermis_triangles else np.zeros(3) for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
topomesh.update_wisp_property('principal_curvature_min',2,np.array([face_principal_curvature_min[t] if t in epidermis_triangles else 0. for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
topomesh.update_wisp_property('principal_curvature_max',2,np.array([face_principal_curvature_max[t] if t in epidermis_triangles else 0. for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
topomesh.update_wisp_property('mean_curvature',2,np.array([(face_principal_curvature_min[t] + face_principal_curvature_max[t])/2. if t in epidermis_triangles else 0. for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
topomesh.update_wisp_property('gaussian_curvature',2,np.array([(face_principal_curvature_min[t]*face_principal_curvature_max[t]) if t in epidermis_triangles else 0. for t in topomesh.wisps(2)]),np.array(list(topomesh.wisps(2))))
elif degree == 3:
if not topomesh.has_wisp_property('principal_curvature_tensor',degree=0,is_computed=True):
compute_topomesh_property(topomesh,'principal_curvature_tensor',degree=0)
if not topomesh.has_wisp_property('epidermis',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',3)
if not topomesh.has_wisp_property('vertices',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'vertices',3)
if not topomesh.has_wisp_property('epidermis',degree=0,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',0)
epidermis_cells = np.array(list(topomesh.wisps(3)))[topomesh.wisp_property('epidermis',3).values(list(topomesh.wisps(3)))]
epidermis_vertices = np.array(list(topomesh.wisps(0)))[topomesh.wisp_property('epidermis',0).values(list(topomesh.wisps(0)))]
cell_epidermis_vertices = np.array([np.intersect1d(topomesh.wisp_property('vertices',degree=3)[c],epidermis_vertices) for c in epidermis_cells])
def row_mean(array):
return np.mean(array,axis=0)
cell_curvature_tensor = array_dict(list(map(row_mean,topomesh.wisp_property('principal_curvature_tensor',degree=0).values(cell_epidermis_vertices))),epidermis_cells)
topomesh.update_wisp_property('principal_curvature_tensor',3,np.array([cell_curvature_tensor[c] if c in epidermis_cells else np.zeros((3,3),float) for c in topomesh.wisps(3)]),np.array(list(topomesh.wisps(3))))
for direction in ['principal_direction_min','principal_direction_max']:
cell_direction = array_dict(list(map(row_mean,topomesh.wisp_property(direction,degree=0).values(cell_epidermis_vertices))),epidermis_cells)
topomesh.update_wisp_property(direction,3,np.array([cell_direction[c] if c in epidermis_cells else np.array([0.,0.,0.]) for c in topomesh.wisps(3)]),np.array(list(topomesh.wisps(3))))
for curvature in ['principal_curvature_min','principal_curvature_max','mean_curvature','gaussian_curvature']:
cell_curvature = array_dict(list(map(np.mean,topomesh.wisp_property(curvature,degree=0).values(cell_epidermis_vertices))),epidermis_cells)
topomesh.update_wisp_property(curvature,3,np.array([cell_curvature[c] if c in epidermis_cells else 0. for c in topomesh.wisps(3)]),np.array(list(topomesh.wisps(3))))
if property_name in ['inertia_tensor','inertia_axes','principal_inertia_axis','inertia_moments']:
assert degree>=2
if degree == 2:
if not topomesh.has_wisp_property('barycenter',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=2)
if not topomesh.has_wisp_property('area',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'area',degree=2)
if not topomesh.has_wisp_property('barycenter',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=1)
if not topomesh.has_wisp_property('length',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'length',degree=1)
if not topomesh.has_wisp_property('edges',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'edges',degree=degree)
face_edge_centers = topomesh.wisp_property('barycenter',1).values(topomesh.wisp_property('edges',2).values(list(topomesh.wisps(2))))
face_edge_weights = topomesh.wisp_property('length',1).values(topomesh.wisp_property('edges',2).values(list(topomesh.wisps(2))))
face_centers = topomesh.wisp_property('barycenter',2).values(list(topomesh.wisps(2)))
face_edge_vectors = np.array([fec - fc for fec,fc in zip(face_edge_centers,face_centers)])
face_inertia_matrix = np.array([np.cov(v,aweights=w,rowvar=False) for v,w in zip(face_edge_vectors,face_edge_weights)])
face_evals, face_evecs = np.linalg.eigh(face_inertia_matrix)
face_inertia_axes = np.array([np.transpose(evec[:,np.argsort(-np.abs(eval))]) for eval,evec in zip(face_evals, face_evecs)])
face_inertia_moments = np.array([eval[np.argsort(-np.abs(eval))] for eval,evec in zip(face_evals, face_evecs)])
topomesh.update_wisp_property('inertia_tensor', 2, dict(zip(topomesh.wisps(2),face_inertia_matrix)))
topomesh.update_wisp_property('inertia_axes', 2, dict(zip(topomesh.wisps(2),face_inertia_axes)))
topomesh.update_wisp_property('principal_inertia_axis', 2, dict(zip(topomesh.wisps(2),face_inertia_axes[:,0])))
topomesh.update_wisp_property('inertia_moments', 2, dict(zip(topomesh.wisps(2),face_inertia_moments)))
if property_name == 'cell_interface':
assert degree==2
if not 'cell_interface' in topomesh.wisp_property_names(degree):
topomesh.add_wisp_property('cell_interface',degree=degree)
if not topomesh.has_wisp_property('cells',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'cells',degree=2)
interface_triangles = np.array(list(topomesh.wisps(2)))[np.where(np.array(list(map(len,topomesh.wisp_property('cells',2).values()))) == 2)[0]]
interface_triangle_cells = np.array([c for c in topomesh.wisp_property('cells',2).values(interface_triangles)])
# def cell_interface(cid1,cid2):
# return topomesh.interface(3,cid1,cid2)
# cell_interfaces = np.array(map(cell_interface,interface_triangle_cells[:,0],interface_triangle_cells[:,1]))
cell_interfaces = vq(np.sort(interface_triangle_cells),np.array(list(topomesh._interface[3].values())))[0]
topomesh.update_wisp_property('cell_interface',degree=2,values=cell_interfaces,keys=interface_triangles)
if property_name == 'interface':
assert degree > 0
if not 'interface' in topomesh.interface_property_names(degree):
topomesh.add_interface_property('interface',degree=degree)
if not topomesh.has_wisp_property('borders',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'borders',degree=degree)
if not topomesh.has_wisp_property('neighbors',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'neighbors',degree=degree)
neighbour_wids = np.array(list(topomesh._interface[degree].values()))
neighbour_borders = topomesh.wisp_property('borders',degree=degree).values(neighbour_wids)
topomesh.update_interface_property('interface',degree=degree,values=np.array(list(map(np.intersect1d,neighbour_borders[:,0],neighbour_borders[:,1]))),keys=np.array(list(topomesh.interfaces(degree))))
if property_name == 'distance':
assert degree > 0
if not 'distance' in topomesh.interface_property_names(degree):
topomesh.add_interface_property('distance',degree=degree)
if not topomesh.has_wisp_property('barycenter',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=degree)
neighbour_wids = np.array(topomesh._interface[degree].values())
neighbour_barycenters = topomesh.wisp_property('barycenter',degree=degree).values(neighbour_wids)
neighbour_vectors = neighbour_barycenters[:,1] - neighbour_barycenters[:,0]
topomesh.update_interface_property('distance',degree=degree,values=np.linalg.norm(neighbour_vectors,axis=1),keys=np.array(list(topomesh.interfaces(degree))))
end_time = time()
if verbose:
print("<-- Computing",property_name,"property (",degree,") [",end_time-start_time,"s]")
[docs]def compute_topomesh_triangle_properties(topomesh,positions=None, verbose=False):
"""Compute an usual set of properties over the faces of a PropertyTopomesh.
The function computes several geometrical properties assuming the structure
passed as argument is a PropertyTopomesh representing a triangular mesh.
The area, perimeter and eccentricity of each triangular face is computed
and the corresponding properties updated in the structure.
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to compute the property.
positions : dict, *optional*
A position dictionary if the property ('barycenter',0) is empty.
Returns
--------
None
Note
--------
The PropertyTopomesh passed as argument is updated.
Example
--------
>>> from cellcomplex.property_topomesh.example_topomesh import square_topomesh
>>> from cellcomplex.property_topomesh.analysis import compute_topomesh_triangle_properties
>>> topomesh = square_topomesh(side_length=1)
>>> compute_topomesh_triangle_properties(topomesh)
>>> print topomesh.wisp_property('area',2)
{0: 0.5, 1: 0.5}
"""
start_time = time()
if verbose: print("--> Computing triangle properties")
if positions is None:
positions = topomesh.wisp_property('barycenter',degree=0)
if not topomesh.has_wisp_property('borders',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'borders',degree=2)
if not topomesh.has_wisp_property('length',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'length',degree=1)
edge_lengths = np.sort(topomesh.wisp_property('length',degree=1).values(topomesh.wisp_property('borders',degree=2).values()))
triangle_perimeters = np.sum(edge_lengths,axis=1)
if not 'perimeter' in topomesh.wisp_property_names(2):
topomesh.add_wisp_property('perimeter',degree=2)
topomesh.update_wisp_property('perimeter',degree=2,values=triangle_perimeters,keys=np.array(list(topomesh.wisps(2))))
triangle_areas = np.sqrt(np.maximum(0,(triangle_perimeters/2.0)*(triangle_perimeters/2.0-edge_lengths[:,0])*(triangle_perimeters/2.0-edge_lengths[:,1])*(triangle_perimeters/2.0-edge_lengths[:,2])))
if not 'area' in topomesh.wisp_property_names(2):
topomesh.add_wisp_property('area',degree=2)
topomesh.update_wisp_property('area',degree=2,values=triangle_areas,keys=np.array(list(topomesh.wisps(2))))
triangle_sinuses = np.zeros_like(edge_lengths,np.float32)
triangle_sinuses[:,0] = np.sqrt(np.maximum(0,np.array(1.0 - np.power(edge_lengths[:,1]**2+edge_lengths[:,2]**2-edge_lengths[:,0]**2,2.0)/np.power(np.maximum(1e-5,2.0*edge_lengths[:,1]*edge_lengths[:,2]),2.0),np.float16)))
triangle_sinuses[:,1] = np.sqrt(np.maximum(0,np.array(1.0 - np.power(edge_lengths[:,2]**2+edge_lengths[:,0]**2-edge_lengths[:,1]**2,2.0)/np.power(np.maximum(1e-5,2.0*edge_lengths[:,2]*edge_lengths[:,0]),2.0),np.float16)))
triangle_sinuses[:,2] = np.sqrt(np.maximum(0,np.array(1.0 - np.power(edge_lengths[:,0]**2+edge_lengths[:,1]**2-edge_lengths[:,2]**2,2.0)/np.power(np.maximum(1e-5,2.0*edge_lengths[:,0]*edge_lengths[:,1]),2.0),np.float16)))
triangle_sinus_eccentricities = 1.0 - (2.0*(triangle_sinuses[:,0]+triangle_sinuses[:,1]+triangle_sinuses[:,2]))/(3*np.sqrt(3))
if not 'eccentricity' in topomesh.wisp_property_names(2):
topomesh.add_wisp_property('eccentricity',degree=2)
#topomesh.update_wisp_property('eccentricity',degree=2,values=triangle_edge_ratios,keys=np.array(list(topomesh.wisps(2))))
#topomesh.update_wisp_property('eccentricity',degree=2,values=triangle_eccentricities,keys=np.array(list(topomesh.wisps(2))))
topomesh.update_wisp_property('eccentricity',degree=2,values=triangle_sinus_eccentricities,keys=np.array(list(topomesh.wisps(2))))
end_time = time()
if verbose: print("<-- Computing triangle properties [",end_time-start_time,"s]")
[docs]def compute_topomesh_vertex_property_from_faces(topomesh, property_name, weighting='area', neighborhood=1, adjacency_sigma=0.5, verbose=False):
"""Compute a property on degree 0 using the same property defined at degree 2.
The vertex property is computed by averaging the properties of its neighbor
faces, weighting them differently according to the chosen method.
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to compute the property.
property_name : str
The name of the property to compute (must be already computed on faces).
weighting : str
The way weights are assigned to each face for the averaging of the property (default is *area*)
* *uniform*: all the faces have the same weight (1)
* *area*: the weight on the faces is equal to their area
* *angle*: the weight of the faces is equal to the incidence angle at the considered vertex (neighborhood=1)
* *cotangent*: the weight of the faces is equal to the sum of cotangent of opposite angles (neighborhood=1)
* *angular sector*: the weight of the faces is computed as the angular sector intersecting face barycenter (neighborhood=1)
neighborhood : int
The ring-distance of faces considered around each vertex (1-ring is immediate neighbor faces).
adjacency_sigma : float
The standard deviation of the gaussian weighting using the ring-distance.
verbose : bool
Whether to display information while the property is being computed
Returns
--------
None
Note
----------
The PropertyTopomesh passed as argument is updated.
Example
----------
>>> from cellcomplex.property_topomesh.example_topomesh import square_topomesh
>>> from cellcomplex.property_topomesh.analysis import compute_topomesh_property
>>> from cellcomplex.property_topomesh.analysis import compute_topomesh_vertex_property_from_faces
>>> topomesh = square_topomesh(side_length=1)
>>> compute_topomesh_property(topomesh,'normal',2,normal_method='orientation')
>>> compute_topomesh_vertex_property_from_faces(topomesh,'normal',weighting='area',neighborhood=3,adjacency_sigma=1.2)
"""
start_time = time()
if verbose: print("--> Computing vertex property from faces")
assert topomesh.has_wisp_property(property_name,2,is_computed=True)
assert weighting in ['uniform','area','angle','cotangent','angular sector']
if neighborhood > 1:
try:
assert weighting in ['uniform','area']
except AssertionError:
raise AssertionError("\""+weighting+"\" weighting can not be used with faces further than 1-ring (neighborhood > 1)")
face_property = topomesh.wisp_property(property_name,2)
assert not face_property.values().dtype == np.object
if neighborhood == 1:
vertex_faces = np.concatenate([list(topomesh.regions(0,v,2)) for v in topomesh.wisps(0)]).astype(np.uint32)
vertex_face_vertices = np.concatenate([v*np.ones_like(list(topomesh.regions(0,v,2))) for v in topomesh.wisps(0)]).astype(np.uint32)
else:
vertex_face_adjacency_matrix = topomesh.nb_wisps(2)*np.ones((max(topomesh.wisps(0))+1,max(topomesh.wisps(2))+1))
vertex_face_adjacency_matrix = topomesh.nb_wisps(2)*np.ones((topomesh.nb_wisps(0)+1,topomesh.nb_wisps(2)+1))
vertex_index = array_dict(dict(zip(topomesh.wisps(0),np.arange(topomesh.nb_wisps(0)))))
face_index = array_dict(dict(zip(topomesh.wisps(2),np.arange(topomesh.nb_wisps(2)))))
vertex_faces = np.concatenate([list(topomesh.regions(0,v,2)) for v in topomesh.wisps(0)]).astype(np.uint32)
vertex_face_vertices = np.concatenate([v*np.ones_like(list(topomesh.regions(0,v,2))) for v in topomesh.wisps(0)]).astype(np.uint32)
compute_topomesh_property(topomesh,'neighbors',2)
vertex_face_adjacency_matrix[(vertex_index.values(vertex_face_vertices),face_index.values(vertex_faces))] = 0
for dist in (np.arange(neighborhood)):
neighbor_face_vertex_indices, neighbor_face_indices = np.where(vertex_face_adjacency_matrix <= dist)
neighbor_face_vertices = np.array(list(topomesh.wisps(0)))[neighbor_face_vertex_indices]
neighbor_faces = np.array(list(topomesh.wisps(2)))[neighbor_face_indices]
neighbor_face_neighbors = np.concatenate(topomesh.wisp_property('neighbors',2).values(neighbor_faces))
neighbor_face_neighbor_vertex = np.concatenate([v*np.ones_like(topomesh.wisp_property('neighbors',2)[f]) for f,v in zip(neighbor_faces,neighbor_face_vertices)])
neighbor_face_neighbor_predecessor = np.concatenate([f*np.ones_like(topomesh.wisp_property('neighbors',2)[f]) for f in neighbor_faces])
current_adjacency = vertex_face_adjacency_matrix[(vertex_index.values(neighbor_face_neighbor_vertex),face_index.values(neighbor_face_neighbors))]
neighbor_adjacency = vertex_face_adjacency_matrix[(vertex_index.values(neighbor_face_neighbor_vertex),face_index.values(neighbor_face_neighbor_predecessor))] + 1
current_adjacency = current_adjacency[np.argsort(-neighbor_adjacency)]
neighbor_face_neighbors = neighbor_face_neighbors [np.argsort(-neighbor_adjacency)]
neighbor_face_neighbor_vertex = neighbor_face_neighbor_vertex[np.argsort(-neighbor_adjacency)]
neighbor_adjacency = neighbor_adjacency[np.argsort(-neighbor_adjacency)]
vertex_face_adjacency_matrix[(vertex_index.values(neighbor_face_neighbor_vertex),face_index.values(neighbor_face_neighbors))] = np.minimum(current_adjacency,neighbor_adjacency)
vertex_face_vertex_indices, vertex_face_indices = np.where(vertex_face_adjacency_matrix < neighborhood)
vertex_face_vertices = np.array(list(topomesh.wisps(0)))[vertex_face_vertex_indices]
vertex_faces = np.array(list(topomesh.wisps(2)))[vertex_face_indices]
vertex_adjacency_gaussian_weight = np.exp(-np.power(vertex_face_adjacency_matrix,2.0)/(2*np.power(adjacency_sigma,2.)))[(vertex_face_vertex_indices, vertex_face_indices)]
if weighting == 'uniform':
vertex_face_weight = np.ones_like(vertex_faces)
elif weighting == 'area':
compute_topomesh_property(topomesh,'area',2)
vertex_face_weight = topomesh.wisp_property('area',2).values(vertex_faces)
elif weighting == 'angle':
compute_topomesh_property(topomesh,'vertices',2)
compute_topomesh_property(topomesh,'angles',2)
vertex_face_face_vertices = topomesh.wisp_property('vertices',2).values(vertex_faces)
vertex_face_face_angles = topomesh.wisp_property('angles',2).values(vertex_faces)
vertex_face_vertex_angles = np.array([angles[vertices == v] for v,vertices,angles in zip(vertex_face_vertices,vertex_face_face_vertices,vertex_face_face_angles)])
if verbose: print(vertex_face_vertex_angles)
vertex_face_weight = vertex_face_vertex_angles[:,0]
elif weighting == 'cotangent':
compute_topomesh_property(topomesh,'vertices',2)
compute_topomesh_property(topomesh,'angles',2)
vertex_face_face_vertices = topomesh.wisp_property('vertices',2).values(vertex_faces)
vertex_face_face_angles = topomesh.wisp_property('angles',2).values(vertex_faces)
vertex_face_opposite_angles = np.array([angles[vertices != v] for v,vertices,angles in zip(vertex_face_vertices,vertex_face_face_vertices,vertex_face_face_angles)])
vertex_face_cotangents = 1./np.tan(vertex_face_opposite_angles)
vertex_face_weight = vertex_face_cotangents.sum(axis=1)
elif weighting == 'angular sector':
compute_topomesh_property(topomesh,'vertices',2)
compute_topomesh_property(topomesh,'angles',2)
compute_topomesh_property(topomesh,'barycenter',2)
vertex_face_face_vertices = topomesh.wisp_property('vertices',2).values(vertex_faces)
vertex_face_face_angles = topomesh.wisp_property('angles',2).values(vertex_faces)
vertex_face_barycenter = topomesh.wisp_property('barycenter',2).values(vertex_faces)
#vertex_positions = topomesh.wisp_property('barycenter',0).values(vertex_face_vertices)
face_vertex_positions = topomesh.wisp_property('barycenter',0).values(vertex_face_face_vertices)
# face_vertex_barycenter_radii = np.linalg.norm(np.array([[vertex_face_barycenter[f] - v for v in face_vertex_positions[f]] for f in vertex_faces]), axis=1)
face_vertex_barycenter_radii = np.linalg.norm(face_vertex_positions - vertex_face_barycenter[:,np.newaxis],axis=1)
vertex_face_weight = np.array(
[r[vertices == v] * angles[vertices == v] for v, vertices, r, angles
in zip(vertex_face_vertices, vertex_face_face_vertices,
face_vertex_barycenter_radii, vertex_face_face_angles)])[:,0]
if neighborhood > 1:
vertex_face_weight = vertex_face_weight*vertex_adjacency_gaussian_weight
vertex_face_property = face_property.values(vertex_faces)
if vertex_face_property.ndim == 1:
vertex_property = nd.sum(vertex_face_weight*vertex_face_property,vertex_face_vertices,index=list(topomesh.wisps(0)))/nd.sum(vertex_face_weight,vertex_face_vertices,index=list(topomesh.wisps(0)))
elif vertex_face_property.ndim == 2:
vertex_property = np.transpose([nd.sum(vertex_face_weight*vertex_face_property[:,k],vertex_face_vertices,index=list(topomesh.wisps(0)))/nd.sum(vertex_face_weight,vertex_face_vertices,index=list(topomesh.wisps(0))) for k in range(vertex_face_property.shape[1])])
elif vertex_face_property.ndim == 3:
vertex_property = np.transpose([[nd.sum(vertex_face_weight*vertex_face_property[:,j,k],vertex_face_vertices,index=list(topomesh.wisps(0)))/nd.sum(vertex_face_weight,vertex_face_vertices,index=list(topomesh.wisps(0))) for k in range(vertex_face_property.shape[2])] for j in range(vertex_face_property.shape[1])])
if property_name in ['normal']:
vertex_property_norm = np.linalg.norm(vertex_property,axis=1)
vertex_property = vertex_property/vertex_property_norm[:,np.newaxis]
vertex_property[np.isnan(vertex_property)] = 0
topomesh.update_wisp_property(property_name,degree=0,values=array_dict(vertex_property,keys=list(topomesh.wisps(0))))
end_time = time()
if verbose: print("<-- Computing vertex property from faces [",end_time-start_time,"s]")
[docs]def compute_topomesh_vertex_property_from_cells(topomesh,property_name,weighting='volume'):
"""Compute a property on degree 0 using the same property defined at degree 3.
The vertex property is computed by averaging the properties of its neighbor
cells, weighting them differently according to the chosen method.
Parameters
--------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to compute the property.
property_name : str
The name of the property to compute (must be already computed on faces).
weighting : str
The way weights are assigned to each cell for the averaging of the property (default is *volume*)
* *uniform*: all the cells have the same weight (1)
* *volume*: the weight on the cells is equal to their volume
Note
--------
The PropertyTopomesh passed as argument is updated.
"""
start_time = time()
print("--> Computing vertex property from cells")
assert topomesh.has_wisp_property(property_name,3,is_computed=True)
assert weighting in ['uniform','volume']
cell_property = topomesh.wisp_property(property_name,3)
assert not cell_property.values().dtype == np.object
vertex_cells = np.concatenate([list(topomesh.regions(0,v,3)) for v in topomesh.wisps(0)]).astype(np.uint32)
vertex_cell_vertices = np.concatenate([v*np.ones_like(list(topomesh.regions(0,v,3))) for v in topomesh.wisps(0)]).astype(np.uint32)
if weighting == 'uniform':
vertex_cell_weight = np.ones_like(vertex_cells)
elif weighting == 'volume':
if not topomesh.has_wisp_property('volume',3,is_computed=True):
compute_topomesh_property(topomesh,'volume',3)
vertex_cell_weight = topomesh.wisp_property('volume',3).values(vertex_cells)
vertex_cell_property = cell_property.values(vertex_cells)
if vertex_cell_property.ndim == 1:
vertex_property = nd.sum(vertex_cell_weight*vertex_cell_property,vertex_cell_vertices,index=list(topomesh.wisps(0)))/nd.sum(vertex_cell_weight,vertex_cell_vertices,index=list(topomesh.wisps(0)))
elif vertex_cell_property.ndim == 2:
vertex_property = np.transpose([nd.sum(vertex_cell_weight*vertex_cell_property[:,k],vertex_cell_vertices,index=list(topomesh.wisps(0)))/nd.sum(vertex_cell_weight,vertex_cell_vertices,index=list(topomesh.wisps(0))) for k in range(vertex_cell_property.shape[1])])
elif vertex_cell_property.ndim == 3:
vertex_property = np.transpose([[nd.sum(vertex_cell_weight*vertex_cell_property[:,j,k],vertex_cell_vertices,index=list(topomesh.wisps(0)))/nd.sum(vertex_cell_weight,vertex_cell_vertices,index=list(topomesh.wisps(0))) for k in range(vertex_cell_property.shape[2])] for j in range(vertex_cell_property.shape[1])])
vertex_property[np.isnan(vertex_property)] = 0
topomesh.update_wisp_property(property_name,degree=0,values=array_dict(vertex_property,keys=list(topomesh.wisps(0))))
end_time = time()
print("<-- Computing vertex property from cells [",end_time-start_time,"s]")
[docs]def compute_topomesh_vertex_property_from_edges(topomesh, property_name, weighting='length', verbose=False):
"""Compute a property on degree 0 using the same property defined at degree 1.
The vertex property is computed by averaging the properties of its neighbor
edges, weighting them differently according to the chosen method.
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to compute the property.
property_name : str
The name of the property to compute (must be already computed on edges).
weighting : str
The way weights are assigned to each edge for the averaging of the property (default is *area*)
* *uniform*: all the edges have the same weight (1)
* *area*: the weight on the edges is equal to their area
* *angle*: the weight of the edges is equal to the incidence angle at the considered vertex (neighborhood=1)
* *cotangent*: the weight of the edges is equal to the sum of cotangent of opposite angles (neighborhood=1)
* *angular sector*: the weight of the edges is computed as the angular sector intersecting edge barycenter (neighborhood=1)
verbose : bool
Whether to display information while the property is being computed
Returns
--------
None
Note
----------
The PropertyTopomesh passed as argument is updated.
"""
start_time = time()
if verbose: print("--> Computing vertex property from edges")
assert topomesh.has_wisp_property(property_name,1,is_computed=True)
assert weighting in ['uniform','length']
edge_property = topomesh.wisp_property(property_name,1)
assert not edge_property.values().dtype == np.object
vertex_edges = np.concatenate([list(topomesh.regions(0,v)) for v in topomesh.wisps(0)]).astype(np.uint32)
vertex_edge_vertices = np.concatenate([v*np.ones_like(list(topomesh.regions(0,v))) for v in topomesh.wisps(0)]).astype(np.uint32)
if weighting == 'uniform':
vertex_edge_weight = np.ones_like(vertex_edges)
elif weighting == 'length':
compute_topomesh_property(topomesh,'length',1)
vertex_edge_weight = topomesh.wisp_property('length',1).values(vertex_edges)
vertex_edge_property = edge_property.values(vertex_edges)
if vertex_edge_property.ndim == 1:
vertex_property = nd.sum(vertex_edge_weight*vertex_edge_property,vertex_edge_vertices,index=list(topomesh.wisps(0)))/nd.sum(vertex_edge_weight,vertex_edge_vertices,index=list(topomesh.wisps(0)))
elif vertex_edge_property.ndim == 2:
vertex_property = np.transpose([nd.sum(vertex_edge_weight*vertex_edge_property[:,k],vertex_edge_vertices,index=list(topomesh.wisps(0)))/nd.sum(vertex_edge_weight,vertex_edge_vertices,index=list(topomesh.wisps(0))) for k in range(vertex_edge_property.shape[1])])
elif vertex_edge_property.ndim == 3:
vertex_property = np.transpose([[nd.sum(vertex_edge_weight*vertex_edge_property[:,j,k],vertex_edge_vertices,index=list(topomesh.wisps(0)))/nd.sum(vertex_edge_weight,vertex_edge_vertices,index=list(topomesh.wisps(0))) for k in range(vertex_edge_property.shape[2])] for j in range(vertex_edge_property.shape[1])])
if property_name in ['normal']:
vertex_property_norm = np.linalg.norm(vertex_property,axis=1)
vertex_property = vertex_property/vertex_property_norm[:,np.newaxis]
vertex_property[np.isnan(vertex_property)] = 0
topomesh.update_wisp_property(property_name,degree=0,values=array_dict(vertex_property,keys=list(topomesh.wisps(0))))
end_time = time()
if verbose: print("<-- Computing vertex property from edges [",end_time-start_time,"s]")
[docs]def compute_topomesh_face_property_from_cells(topomesh,property_name,weighting='uniform'):
"""Compute a property on degree 2 using the same property defined at degree 3.
The face property is computed by averaging the properties of its neighbor
cells, weighting them differently according to the chosen method.
Parameters
--------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to compute the property.
property_name : str
The name of the property to compute (must be already computed on faces).
weighting : str
The way weights are assigned to each cell for the averaging of the property (default is *volume*)
* *uniform*: all the cells have the same weight (1)
* *volume*: the weight on the cells is equal to their volume
Returns
--------
None
Note
--------
The PropertyTopomesh passed as argument is updated.
"""
start_time = time()
print("--> Computing face property from cells")
assert topomesh.has_wisp_property(property_name,3,is_computed=True)
assert weighting in ['uniform','volume']
cell_property = topomesh.wisp_property(property_name,3)
assert not cell_property.values().dtype == np.object
face_cells = np.concatenate([list(topomesh.regions(2,f)) for f in topomesh.wisps(2)]).astype(np.uint32)
face_cell_faces = np.concatenate([f*np.ones_like(list(topomesh.regions(2,f))) for f in topomesh.wisps(2)]).astype(np.uint32)
if weighting == 'uniform':
face_cell_weight = np.ones_like(face_cells)
elif weighting == 'volume':
if not topomesh.has_wisp_property('volume',3,is_computed=True):
compute_topomesh_property(topomesh,'volume',3)
face_cell_weight = topomesh.wisp_property('volume',3).values(face_cells)
face_cell_property = cell_property.values(face_cells)
if face_cell_property.ndim == 1:
face_property = nd.sum(face_cell_weight*face_cell_property,face_cell_faces,index=list(topomesh.wisps(2)))/nd.sum(face_cell_weight,face_cell_faces,index=list(topomesh.wisps(2)))
elif face_cell_property.ndim == 2:
face_property = np.transpose([nd.sum(face_cell_weight*face_cell_property[:,k],face_cell_faces,index=list(topomesh.wisps(2)))/nd.sum(face_cell_weight,face_cell_faces,index=list(topomesh.wisps(2))) for k in range(face_cell_property.shape[1])])
elif face_cell_property.ndim == 3:
face_property = np.transpose([[nd.sum(face_cell_weight*face_cell_property[:,j,k],face_cell_faces,index=list(topomesh.wisps(2)))/nd.sum(face_cell_weight,face_cell_faces,index=list(topomesh.wisps(2))) for k in range(face_cell_property.shape[2])] for j in range(face_cell_property.shape[1])])
face_property[np.isnan(face_property)] = 0
topomesh.update_wisp_property(property_name,degree=2,values=array_dict(face_property,keys=list(topomesh.wisps(2))))
end_time = time()
print("<-- Computing face property from cells [",end_time-start_time,"s]")
[docs]def compute_topomesh_cell_property_from_faces(topomesh, property_name, reduce='mean', weighting='area', verbose=False):
"""Compute a property on degree 3 using the same property defined at degree 2.
The cell property is computed by averaging or summing the properties of its
border faces, weighting them differently according to the chosen method.
Parameters
--------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to compute the property.
property_name : str
The name of the property to compute (must be already computed on faces).
reduce : str
The way weighted face properties are combined to compute the cell property (default is *mean*):
* *mean*: the cell property is the weighted average of face properties
* *sum*: the cell property is the weighted sum of face properties
weighting : str
The way weights are assigned to each face to compute the property (default is *area*):
* *uniform*: all the faces have the same weight (1)
* *area*: the weight on the faces is equal to their area
Returns
--------
None
Note
--------
The PropertyTopomesh passed as argument is updated.
"""
start_time = time()
if verbose: print("--> Computing cell property from faces")
assert topomesh.has_wisp_property(property_name,2,is_computed=True)
assert weighting in ['uniform','area']
face_property = topomesh.wisp_property(property_name,2)
cell_faces = np.concatenate([list(topomesh.borders(3,c)) for c in topomesh.wisps(3)]).astype(np.uint32)
cell_face_cells = np.concatenate([c*np.ones(topomesh.nb_borders(3,c)) for c in topomesh.wisps(3)]).astype(np.uint32)
if weighting == 'uniform':
cell_face_weight = np.ones_like(cell_faces)
elif weighting == 'area':
compute_topomesh_property(topomesh,'area',2)
cell_face_weight = topomesh.wisp_property('area',2).values(cell_faces)
cell_face_property = face_property.values(cell_faces)
if cell_face_property.ndim == 1:
if reduce in ['mean', 'sum']:
cell_property = nd.sum(cell_face_weight*cell_face_property,cell_face_cells,index=list(topomesh.wisps(3)))
if reduce == 'mean':
cell_property = cell_property/nd.sum(cell_face_weight,cell_face_cells,index=list(topomesh.wisps(3)))
elif cell_face_property.ndim == 2:
if reduce in ['mean', 'sum']:
cell_property = np.transpose([nd.sum(cell_face_weight*cell_face_property[:,k],cell_face_cells,index=list(topomesh.wisps(3))) for k in range(cell_face_property.shape[1])])
if reduce == 'mean':
cell_property = cell_property / nd.sum(cell_face_weight,cell_face_cells,index=list(topomesh.wisps(3)))[:,np.newaxis]
elif cell_face_property.ndim == 3:
if reduce in ['mean', 'sum']:
cell_property = np.transpose([[nd.sum(cell_face_weight*cell_face_property[:,j,k],cell_face_cells,index=list(topomesh.wisps(3))) for k in range(cell_face_property.shape[2])] for j in range(cell_face_property.shape[1])])
if reduce == 'mean':
cell_property = cell_property / nd.sum(cell_face_weight,cell_face_cells,index=list(topomesh.wisps(3)))[:,np.newaxis,np.newaxis]
topomesh.update_wisp_property(property_name,degree=3,values=array_dict(cell_property,keys=list(topomesh.wisps(3))))
end_time = time()
if verbose: print("<-- Computing cell property from faces [",end_time-start_time,"s]")
[docs]def topomesh_property_gaussian_filtering(topomesh,property_name,degree,neighborhood=3,adjacency_sigma=1.0,distance_sigma=1.0):
"""Filter an existing property by a gaussian-like operator.
A new value of an existing property is computed for each element of a given
degree as a linear combination of all the property values. The coefficient
associated with each other element is defined relatively to its distance
in both space and topology (ring-distance) using the product of 2 gaussian
functions.
Parameters
--------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to compute the property.
property_name : str
The name of the property to compute (must be already computed on faces).
degree : int
The degree of the elements on which to compute the property.
neighborhood : int
The maximal ring-distance up to which elements have a non-zero ceofficient.
adjacency_sigma : float
The standard deviation of the gaussian function based on the ring-distance.
distance_sigma : float
The standard deviation of the gaussian function based on the Euclidean distance.
Returns
--------
None
Note
--------
The PropertyTopomesh passed as argument is updated.
"""
assert topomesh.has_wisp_property(property_name,degree=degree,is_computed=True)
if not topomesh.has_wisp_property('barycenter',degree=degree,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=degree)
adjacency_matrix = topomesh.nb_wisps(degree)*np.ones((topomesh.nb_wisps(degree),topomesh.nb_wisps(degree)),int)
adjacency_matrix[tuple([np.arange(topomesh.nb_wisps(degree)),np.arange(topomesh.nb_wisps(degree))])] = 0
wisp_index = array_dict(np.arange(topomesh.nb_wisps(degree)),list(topomesh.wisps(degree)))
for i,wid in enumerate(topomesh.wisps(degree)):
if degree>1:
adjacency_matrix[tuple([i*np.ones(topomesh.nb_border_neighbors(degree,wid),int),wisp_index.values(list(topomesh.border_neighbors(degree,wid)))])] = 1
else:
adjacency_matrix[tuple([i*np.ones(topomesh.nb_region_neighbors(degree,wid),int),wisp_index.values(list(topomesh.region_neighbors(degree,wid)))])] = 1
for dist in (np.arange(neighborhood)+1):
for i,wid in enumerate(topomesh.wisps(degree)):
neighbors = np.where(adjacency_matrix[i] <= dist)[0]
neighbor_neighbors = np.where(adjacency_matrix[neighbors] <= dist)
neighbor_neighbor_successor = neighbor_neighbors[1]
neighbor_neighbor_predecessor = neighbors[neighbor_neighbors[0]]
current_adjacency = adjacency_matrix[i][neighbor_neighbor_successor]
neighbor_adjacency = adjacency_matrix[i][neighbor_neighbor_predecessor] + adjacency_matrix[neighbors][neighbor_neighbors]
current_adjacency = current_adjacency[np.argsort(-neighbor_adjacency)]
neighbor_neighbor_successor = neighbor_neighbor_successor[np.argsort(-neighbor_adjacency)]
neighbor_adjacency = neighbor_adjacency[np.argsort(-neighbor_adjacency)]
adjacency_matrix[i][neighbor_neighbor_successor] = np.minimum(current_adjacency,neighbor_adjacency)
distance_matrix = np.array([vq(topomesh.wisp_property('barycenter',degree).values(),np.array([topomesh.wisp_property('barycenter',degree)[v]]))[1] for v in topomesh.wisps(degree)])
vertex_adjacency_gaussian = np.exp(-np.power(adjacency_matrix,2.0)/(2*np.power(adjacency_sigma,2.)))
vertex_distance_gaussian = np.exp(-np.power(distance_matrix,2.0)/(2*np.power(distance_sigma,2.)))
vertex_gaussian = vertex_adjacency_gaussian*vertex_distance_gaussian
properties = topomesh.wisp_property(property_name,degree=degree).values(list(topomesh.wisps(degree)))
if properties.ndim == 1:
filtered_properties = (vertex_gaussian*properties[:,np.newaxis]).sum(axis=0) / vertex_gaussian.sum(axis=0)
else:
filtered_properties = (vertex_gaussian[...,np.newaxis]*properties[:,np.newaxis]).sum(axis=0) / vertex_gaussian.sum(axis=0)[...,np.newaxis]
topomesh.update_wisp_property(property_name,degree=degree,values=filtered_properties,keys=np.array(list(topomesh.wisps(degree))))
[docs]def is_triangular(topomesh):
"""Check wether the structure is a triangular mesh.
Look at the number of vertices of each face of the topomesh to decide
whether it can be processed as a triangular mesh in all the algorithms.
Parameters
--------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure to evaluate.
Returns
--------
triangular : bool
True if the mesh has only triangular faces, False otherwise
"""
if topomesh.nb_wisps(2)==0:
return True
else:
compute_topomesh_property(topomesh,'vertices',2)
face_vertices = topomesh.wisp_property('vertices',2).values()
return (face_vertices.ndim == 2) and (face_vertices.shape[1]==3)
