# -*- coding: utf-8 -*-
# -*- python -*-
#
# PropertyTopomesh
#
# Copyright 2015-2016 INRIA - CIRAD - INRA
#
# File author(s): Guillaume Cerutti <guillaume.cerutti@inria.fr>
#
# File contributor(s): Guillaume Cerutti <guillaume.cerutti@inria.fr>
#
# 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 Optimization
=============================
.. autosummary::
:toctree: generated/
property_topomesh_vertices_deformation -- Geometrical multi-criterion optimization
property_topomesh_edge_flip_optimization -- Topological optimization by successive edge flips
property_topomesh_edge_split_optimization -- Topological optimization by successive edge splits
property_topomesh_edge_collapse_optimization -- Topological optimization by successive edge collapses
property_topomesh_isotropic_remeshing -- Global triangle quality optimization
topomesh_triangle_split -- Simple remeshing by global subdivision of triangles
"""
from time import time
from copy import deepcopy
import logging
import numpy as np
from scipy import ndimage as nd
from scipy.cluster.vq import vq
from cellcomplex.property_topomesh import PropertyTopomesh
from cellcomplex.property_topomesh.analysis import is_triangular, compute_topomesh_property, compute_topomesh_triangle_properties
from cellcomplex.property_topomesh.topological_operations import topomesh_flip_edge, topomesh_split_edge, topomesh_collapse_edge
from cellcomplex.utils.array_dict import array_dict
from cellcomplex.property_topomesh.utils.geometry_tools import triangle_geometric_features
def _property_topomesh_cell_vertex_force(topomesh,gradient_derivatives,voxelsize, verbose=False, debug=False, loglevel=0):
"""
Compute for each vertex of the topomesh the force guiding its displacement towards a cell vertex
"""
gradient_x = gradient_derivatives[0]
gradient_y = gradient_derivatives[1]
gradient_z = gradient_derivatives[2]
vertices_coords = np.rollaxis(np.array(topomesh.wisp_property('barycenter',degree=0).values()/np.array(voxelsize),np.uint16),1)
vertices_coords = np.maximum(np.minimum(vertices_coords,([gradient_x.shape[0]-1],[gradient_x.shape[1]-1],[gradient_x.shape[2]-1])),0)
# logging.info("".join([" " for l in range(loglevel)])+vertices_coords
gradient_force = np.rollaxis(np.array([gradient_x[tuple(vertices_coords)],gradient_y[tuple(vertices_coords)],gradient_z[tuple(vertices_coords)]]),1)*np.array(voxelsize)[np.newaxis,:]
return gradient_force
def _property_topomesh_triangle_regularization_force(topomesh, verbose=False, debug=False, loglevel=0):
"""todo"""
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('length',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'length',degree=1)
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values()
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.concatenate([triangle_vertices,rotated_triangle_vertices,antirotated_triangle_vertices],axis=0)
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.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_areas[np.where(triangle_areas==0.0)] = 0.001
average_area = np.mean(triangle_areas)
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
sinus_edge_1_force = (((np.power(triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2,2.0)-np.power(triangle_edge_lengths[:,1],4.0))*triangle_edge_lengths[:,2])/(2.0*np.power(triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2],3.0)))/(2.0*triangle_sinuses[:,0])
sinus_edge_1_force += ((2.*(triangle_edge_lengths[:,2]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,1]**2)*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2])/(2.0*np.power(triangle_edge_lengths[:,0]*triangle_edge_lengths[:,2],3.0)))/(2.0*triangle_sinuses[:,1])
sinus_edge_1_force += (((np.power(triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2,2.0)-np.power(triangle_edge_lengths[:,1],4.0))*triangle_edge_lengths[:,0])/(2.0*np.power(triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1],3.0)))/(2.0*triangle_sinuses[:,2])
sinus_edge_2_force = (((np.power(triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,0]**2,2.0)-np.power(triangle_edge_lengths[:,2],4.0))*triangle_edge_lengths[:,1])/(2.0*np.power(triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2],3.0)))/(2.0*triangle_sinuses[:,0])
sinus_edge_2_force += (((np.power(triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,0]**2,2.0)-np.power(triangle_edge_lengths[:,2],4.0))*triangle_edge_lengths[:,0])/(2.0*np.power(triangle_edge_lengths[:,0]*triangle_edge_lengths[:,2],3.0)))/(2.0*triangle_sinuses[:,1])
sinus_edge_2_force += ((2.*(triangle_edge_lengths[:,1]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,2]**2)*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2])/(2.0*np.power(triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1],3.0)))/(2.0*triangle_sinuses[:,2])
#sinus_unitary_force = -(sinus_edge_1_force[:,np.newaxis]*triangle_edge_directions[:,1] + sinus_edge_2_force[:,np.newaxis]*triangle_edge_directions[:,2])*triangle_areas[:,np.newaxis]
sinus_unitary_force = -(sinus_edge_1_force[:,np.newaxis]*triangle_edge_vectors[:,1] + sinus_edge_2_force[:,np.newaxis]*triangle_edge_vectors[:,2])
area_edge_1_force = -((triangle_areas - average_area)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,1]**2)*triangle_edge_lengths[:,1])/(4.0*triangle_areas)
area_edge_2_force = -((triangle_areas - average_area)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)*triangle_edge_lengths[:,2])/(4.0*triangle_areas)
#area_edge_1_force = -((triangle_areas - average_area)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,1]**2)*triangle_edge_lengths[:,1])/(4.0*average_area)
#area_edge_2_force = -((triangle_areas - average_area)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)*triangle_edge_lengths[:,2])/(4.0*average_area)
area_unitary_force = -(area_edge_1_force[:,np.newaxis]*triangle_edge_directions[:,1] + area_edge_2_force[:,np.newaxis]*triangle_edge_directions[:,2])
#triangle_unitary_force = 8.0*sinus_unitary_force
#triangle_unitary_force = sinus_unitary_force
triangle_unitary_force = sinus_unitary_force*np.power(average_area,1)
#triangle_unitary_force = sinus_unitary_force*np.power(average_area,1) + area_unitary_force
#triangle_unitary_force = sinus_unitary_force + area_unitary_force
#triangle_unitary_force = 0.02*area_unitary_force + 8.0*sinus_unitary_force
# triangle_unitary_force = 0.02*area_unitary_force
triangle_force = np.transpose([nd.sum(triangle_unitary_force[:,k],triangle_vertices[:,0],index=list(topomesh.wisps(0))) for k in range(3)])
return triangle_force
def _property_topomesh_area_smoothing_force(topomesh,target_areas=None, verbose=False, debug=False, loglevel=0):
compute_topomesh_property(topomesh,'vertices',degree=2)
compute_topomesh_property(topomesh,'length',degree=1)
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values()
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_areas[np.where(triangle_areas==0.0)] = 0.001
if target_areas is None:
target_areas = triangle_areas.mean()*np.ones(len(triangle_areas))
else:
target_areas = np.tile(target_areas,(3))
area_edge_1_force = (target_areas-triangle_areas)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,1]**2)*triangle_edge_lengths[:,1]
area_edge_2_force = (target_areas-triangle_areas)*(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)*triangle_edge_lengths[:,2]
area_unitary_force = -(area_edge_1_force[:,np.newaxis]*triangle_edge_directions[:,1] + area_edge_2_force[:,np.newaxis]*triangle_edge_directions[:,2])
triangle_force = np.transpose([nd.sum(area_unitary_force[:,d],triangle_vertices[:,0],index=list(topomesh.wisps(0))) for d in range(3)])
return triangle_force
def _property_topomesh_laplacian_smoothing_force(topomesh,cellwise_smoothing=False, verbose=False, debug=False, loglevel=0):
if not topomesh.has_wisp_property('vertices',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=1)
edge_vertices = topomesh.wisp_property('vertices',degree=1).values()
reversed_edge_vertices = np.transpose([edge_vertices[:,1],edge_vertices[:,0]])
edge_vertices = np.append(edge_vertices,reversed_edge_vertices,axis=0)
if not topomesh.has_wisp_property('valency',degree=0,is_computed=True):
compute_topomesh_property(topomesh,'valency',degree=0)
vertices_degrees = topomesh.wisp_property('valency',degree=0).values()
if cellwise_smoothing:
laplacian_force = np.zeros_like(topomesh.wisp_property('barycenter',degree=0).values(),np.float32)
for c in topomesh.wisps(3):
if not topomesh.has_wisp_property('vertices',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=3)
cell_vertices = topomesh.wisp_property('vertices',degree=3)[c]
cell_edges = np.sum(nd.sum(np.ones_like(cell_vertices),cell_vertices,index=edge_vertices),axis=1)
cell_edge_vertices = edge_vertices[np.where(cell_edges==2)[0]]
cell_edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(cell_edge_vertices[:,1]) - topomesh.wisp_property('barycenter',degree=0).values(cell_edge_vertices[:,0])
cell_laplacian_force = np.transpose([nd.sum(cell_edge_vectors[:,0],cell_edge_vertices[:,0],index=cell_vertices),
nd.sum(cell_edge_vectors[:,1],cell_edge_vertices[:,0],index=cell_vertices),
nd.sum(cell_edge_vectors[:,2],cell_edge_vertices[:,0],index=cell_vertices)])
laplacian_force[cell_vertices] += cell_laplacian_force/vertices_degrees[cell_vertices,np.newaxis]
else:
edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,1]) - topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,0])
laplacian_force = np.transpose([nd.sum(edge_vectors[:,d],edge_vertices[:,0],index=list(topomesh.wisps(0))) for d in [0,1,2]])
laplacian_force = laplacian_force/vertices_degrees[:,np.newaxis]
return laplacian_force
def _property_topomesh_cotangent_laplacian_smoothing_force(topomesh, verbose=False, debug=False, loglevel=0):
compute_topomesh_property(topomesh,'vertices',degree=2)
compute_topomesh_property(topomesh,'length',degree=1)
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, 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.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_cosines[np.where(np.abs(triangle_cosines) < np.power(10.,-6))] = np.power(10.,-6)
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 < np.power(10.,-6))] = np.power(10.,-6)
triangle_cotangent_weights = (triangle_cosines/triangle_sinuses)
#triangle_cotangent_weights = 1./triangle_edge_lengths
#triangle_cotangent_weights = 0.5*(triangle_cosines/triangle_sinuses)/(triangle_edge_lengths + np.power(10.,-10))
#triangle_cotangent_vectors = (triangle_cosines/triangle_sinuses)[...,np.newaxis] * triangle_edge_vectors/triangle_edge_lengths[:,np.newaxis]
triangle_cotangent_vectors = triangle_cotangent_weights[...,np.newaxis] * triangle_edge_vectors
#triangle_cotangent_vectors = triangle_cotangent_weights[...,np.newaxis] * triangle_edge_directions
# triangle_cotangent_vectors = 1./np.tan(triangle_angles)[...,np.newaxis] * triangle_edge_vectors
vertex_cotangent_force = np.transpose([nd.sum(triangle_cotangent_vectors[:,1,d]+triangle_cotangent_vectors[:,2,d],triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0)))) for d in range(3)])
vertex_cotangent_sum = nd.sum(triangle_cotangent_weights[:,1] + triangle_cotangent_weights[:,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))))
#laplacian_force = vertex_cotangent_force
laplacian_force = vertex_cotangent_force/vertex_cotangent_sum[...,np.newaxis]
#laplacian_force = vertex_cotangent_force/(4.*vertex_area[...,np.newaxis])
laplacian_force[np.where(np.isnan(laplacian_force))] = 0.0
laplacian_force[np.where(np.isinf(laplacian_force))] = 0.0
return laplacian_force
def _property_topomesh_gaussian_smoothing_force(topomesh,gaussian_sigma=10.0, verbose=False, debug=False, loglevel=0):
compute_topomesh_property(topomesh,'vertices',degree=1)
edge_vertices = topomesh.wisp_property('vertices',degree=1).values()
reversed_edge_vertices = np.transpose([edge_vertices[:,1],edge_vertices[:,0]])
edge_vertices = np.append(edge_vertices,reversed_edge_vertices,axis=0)
edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,1]) - topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,0])
edge_lengths = np.linalg.norm(edge_vectors,axis=1)
gaussian_edge_lengths = np.exp(-np.power(edge_lengths,2.0)/np.power(gaussian_sigma,2.0))
gaussian_edge_vectors = gaussian_edge_lengths[:,np.newaxis]*edge_vectors
gaussian_force = np.transpose([nd.sum(gaussian_edge_vectors[:,d],edge_vertices[:,0],index=list(topomesh.wisps(0))) for d in [0,1,2]])
vertices_weights = 1.+nd.sum(gaussian_edge_lengths,edge_vertices[:,0],index=list(topomesh.wisps(0)))
gaussian_force = gaussian_force/vertices_weights[:,np.newaxis]
return gaussian_force
def _property_topomesh_taubin_smoothing_force(topomesh,gaussian_sigma=10.0,positive_factor=0.33,negative_factor=-0.34,cellwise_smoothing=True, verbose=False, debug=False, loglevel=0):
if not topomesh.has_wisp_property('vertices',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'vertices',degree=1)
if cellwise_smoothing:
edge_vertices = topomesh.wisp_property('vertices',degree=1).values(np.concatenate([np.array(list(topomesh.borders(3,c,2)),int) for c in topomesh.wisps(3)]))
else:
edge_vertices = topomesh.wisp_property('vertices',degree=1).values()
reversed_edge_vertices = np.transpose([edge_vertices[:,1],edge_vertices[:,0]])
edge_vertices = np.append(edge_vertices,reversed_edge_vertices,axis=0)
edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,1]) - topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,0])
edge_lengths = np.linalg.norm(edge_vectors,axis=1)
gaussian_edge_lengths = np.exp(-np.power(edge_lengths,2.0)/np.power(gaussian_sigma,2.0))
gaussian_edge_vectors = gaussian_edge_lengths[:,np.newaxis]*edge_vectors
taubin_force = np.transpose([nd.sum(gaussian_edge_vectors[:,0],edge_vertices[:,0],index=list(topomesh.wisps(0))),
nd.sum(gaussian_edge_vectors[:,1],edge_vertices[:,0],index=list(topomesh.wisps(0))),
nd.sum(gaussian_edge_vectors[:,2],edge_vertices[:,0],index=list(topomesh.wisps(0)))])
#vertices_weights = 1.+nd.sum(gaussian_edge_lengths,edge_vertices[:,0],index=list(topomesh.wisps(0)))
vertices_weights = nd.sum(gaussian_edge_lengths,edge_vertices[:,0],index=list(topomesh.wisps(0)))
taubin_positive_force = positive_factor*taubin_force/vertices_weights[:,np.newaxis]
taubin_positions = array_dict(topomesh.wisp_property('barycenter',degree=0).values(list(topomesh.wisps(0)))+taubin_positive_force, list(topomesh.wisps(0)))
edge_vectors = taubin_positions.values(edge_vertices[:,1]) - taubin_positions.values(edge_vertices[:,0])
edge_lengths = np.linalg.norm(edge_vectors,axis=1)
gaussian_edge_lengths = np.exp(-np.power(edge_lengths,2.0)/np.power(gaussian_sigma,2.0))
gaussian_edge_vectors = gaussian_edge_lengths[:,np.newaxis]*edge_vectors
taubin_force = np.transpose([nd.sum(gaussian_edge_vectors[:,0],edge_vertices[:,0],index=list(topomesh.wisps(0))),
nd.sum(gaussian_edge_vectors[:,1],edge_vertices[:,0],index=list(topomesh.wisps(0))),
nd.sum(gaussian_edge_vectors[:,2],edge_vertices[:,0],index=list(topomesh.wisps(0)))])
#vertices_weights = 1.+nd.sum(gaussian_edge_lengths,edge_vertices[:,0],index=list(topomesh.wisps(0)))
vertices_weights = nd.sum(gaussian_edge_lengths,edge_vertices[:,0],index=list(topomesh.wisps(0)))
taubin_negative_force = negative_factor*taubin_force/vertices_weights[:,np.newaxis]
taubin_force = taubin_positive_force + taubin_negative_force
return taubin_force
def _property_topomesh_mean_curvature_smoothing_force(topomesh, verbose=False, debug=False, loglevel=0):
"""todo"""
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)
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, 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.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])
return vertex_mean_curvature_vectors
def _property_topomesh_laplacian_epidermis_convexity_force(topomesh, verbose=False, debug=False, loglevel=0):
"""todo"""
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('barycenter',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=3)
if not topomesh.has_wisp_property('epidermis',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',degree=1)
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('cells',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'cells',degree=1)
epidermis_convexity_force = array_dict(np.zeros_like(topomesh.wisp_property('barycenter',degree=0).values(),np.float32),np.array(list(topomesh.wisps(0))))
edge_cells_degree = np.array(list(map(len,topomesh.wisp_property('cells',degree=1).values())))
# edge_vertices = topomesh.wisp_property('vertices',degree=1).values()[np.where((topomesh.wisp_property('epidermis',degree=1).values()))]
# edge_vertices = topomesh.wisp_property('vertices',degree=1).values()[np.where((topomesh.wisp_property('epidermis',degree=1).values())&(edge_cells_degree>1))]
edge_vertices = topomesh.wisp_property('vertices',degree=1).values()[np.where((edge_cells_degree>2)|((topomesh.wisp_property('epidermis',degree=1).values())&(edge_cells_degree>1)))]
epidermis_vertices = np.unique(edge_vertices)
vertices_degrees = np.array(nd.sum(np.ones_like(edge_vertices),edge_vertices,index=epidermis_vertices),int)
reversed_edge_vertices = np.transpose([edge_vertices[:,1],edge_vertices[:,0]])
edge_vertices = np.append(edge_vertices,reversed_edge_vertices,axis=0)
edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,1]) - topomesh.wisp_property('barycenter',degree=0).values(edge_vertices[:,0])
laplacian_force = np.transpose([nd.sum(edge_vectors[:,0],edge_vertices[:,0],index=epidermis_vertices),
nd.sum(edge_vectors[:,1],edge_vertices[:,0],index=epidermis_vertices),
nd.sum(edge_vectors[:,2],edge_vertices[:,0],index=epidermis_vertices)])
laplacian_force[vertices_degrees>0] = laplacian_force[vertices_degrees>0]/vertices_degrees[vertices_degrees>0][:,np.newaxis]
# vertex_cell_barycenter = np.array([np.mean(topomesh.wisp_property('barycenter',degree=3).values(c),axis=0) for c in topomesh.wisp_property('cells',0).values(epidermis_vertices)])
# vertices_directions = topomesh.wisp_property('barycenter',degree=0).values(epidermis_vertices) - vertex_cell_barycenter
# vertices_directions = vertices_directions/np.linalg.norm(vertices_directions,axis=1)[:,np.newaxis]
# vertices_products = np.einsum('ij,ij->i',laplacian_force,vertices_directions)
# convex_points = np.where(vertices_products<0.)[0]
# laplacian_force[convex_points] -= vertices_directions[convex_points]*vertices_products[convex_points,np.newaxis]
# laplacian_force[convex_points] = np.array([0.,0.,0.])
epidermis_convexity_force.update(laplacian_force,keys=epidermis_vertices,erase_missing_keys=False)
return epidermis_convexity_force.values()
def _property_topomesh_interface_planarization_force(topomesh, verbose=False, debug=False, loglevel=0):
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('vertices',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'vertices',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('cell_interface',degree=2,is_computed=False):
compute_topomesh_property(topomesh,'cell_interface',degree=2)
if not topomesh.has_interface_property('interface',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'interface',degree=3)
interface_planarization_force = np.zeros_like(topomesh.wisp_property('barycenter',degree=0).values(),np.float32)
cell_interfaces = topomesh.interface_property('interface',degree=3).keys()
interface_triangles = topomesh.wisp_property('cell_interface',degree=2).keys()
triangle_cell_interface = topomesh.wisp_property('cell_interface',degree=2).values()
triangle_cells = np.array(array_dict(topomesh._interface[3]).values(triangle_cell_interface))
triangle_cell_barycenters = topomesh.wisp_property('barycenter',degree=3).values(triangle_cells)
triangle_interface_directions = triangle_cell_barycenters[:,1] - triangle_cell_barycenters[:,0]
triangle_normal_vectors = topomesh.wisp_property('normal',degree=2).values(interface_triangles)
# reversed_normals = np.where(np.dot(triangle_normal_vectors,triangle_interface_directions) < 0)[0]
reversed_normals = np.where(np.einsum('ij,ij->i',triangle_normal_vectors,triangle_interface_directions) < 0)[0]
triangle_normal_vectors[reversed_normals] = -triangle_normal_vectors[reversed_normals]
interface_normal_vectors = np.transpose([nd.mean(triangle_normal_vectors[:,0],triangle_cell_interface,index=cell_interfaces),
nd.mean(triangle_normal_vectors[:,1],triangle_cell_interface,index=cell_interfaces),
nd.mean(triangle_normal_vectors[:,2],triangle_cell_interface,index=cell_interfaces)])
interface_normal_vectors = interface_normal_vectors/np.linalg.norm(interface_normal_vectors,axis=1)[:,np.newaxis]
triangle_target_normals = interface_normal_vectors[triangle_cell_interface]
triangle_barycenters = topomesh.wisp_property('barycenter',degree=2).values(interface_triangles)
interface_barycenters = np.transpose([nd.mean(triangle_barycenters[:,0],triangle_cell_interface,index=cell_interfaces),
nd.mean(triangle_barycenters[:,1],triangle_cell_interface,index=cell_interfaces),
nd.mean(triangle_barycenters[:,2],triangle_cell_interface,index=cell_interfaces)])
triangle_interface_barycenters = interface_barycenters[triangle_cell_interface]
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(interface_triangles)
interface_vertices = np.unique(triangle_vertices)
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_point_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_vertices[:,0]) - np.tile(triangle_interface_barycenters,(3,1))
# triangle_point_projectors = -np.dot(triangle_point_vectors,triange_target_normals)[:,np.newaxis]*triangle_target_normals
# logging.info("".join([" " for l in range(loglevel)])+triangle_point_vectors.shape,triangle_target_normals.shape
triangle_point_projectors = -np.einsum('ij,ij->i',triangle_point_vectors,np.tile(triangle_target_normals,(3,1)))[:,np.newaxis]*np.tile(triangle_target_normals,(3,1))
# logging.info("".join([" " for l in range(loglevel)])+triangle_point_projectors.shape
triangle_projectors_norm = np.linalg.norm(triangle_point_projectors,axis=1)
non_planar_points = np.where(triangle_projectors_norm > 1.)[0]
triangle_point_projectors[non_planar_points] = triangle_point_projectors[non_planar_points]/triangle_projectors_norm[non_planar_points,np.newaxis]
# planarization_unitary_force = normal_errors[:,np.newaxis]*point_projectors
planarization_unitary_force = triangle_point_projectors
interface_planarization_force = np.transpose([nd.sum(planarization_unitary_force[:,0],triangle_vertices[:,0],index=list(topomesh.wisps(0))),
nd.sum(planarization_unitary_force[:,1],triangle_vertices[:,0],index=list(topomesh.wisps(0))),
nd.sum(planarization_unitary_force[:,2],triangle_vertices[:,0],index=list(topomesh.wisps(0)))])
interface_planarization_force[np.isnan(interface_planarization_force)] = 0.
return interface_planarization_force
def _property_topomesh_epidermis_convexity_force(topomesh, verbose=False, debug=False, loglevel=0):
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('vertices',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'vertices',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('cells',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'cells',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)))[np.where(topomesh.wisp_property('epidermis',2).values())]
triangle_cell = np.array([c[0] for c in topomesh.wisp_property('cells',2).values(epidermis_triangles)])
epidermis_cells = np.unique(triangle_cell)
triangle_cell_barycenters = topomesh.wisp_property('barycenter',degree=3).values(triangle_cell)
triangle_barycenters = topomesh.wisp_property('barycenter',degree=2).values(epidermis_triangles)
cell_epidermis_barycenters = np.transpose([nd.sum(triangle_barycenters[:,k],triangle_cell,index=list(topomesh.wisps(3))) for k in range(3)])
cell_epidermis_n_triangles = nd.sum(np.ones_like(triangle_cell),triangle_cell,index=list(topomesh.wisps(3)))
cell_epidermis_barycenters[cell_epidermis_n_triangles>0] /= cell_epidermis_n_triangles[cell_epidermis_n_triangles>0][:,np.newaxis]
cell_epidermis_barycenters = array_dict(cell_epidermis_barycenters,keys=list(topomesh.wisps(3)))
triangle_epidermis_barycenters = cell_epidermis_barycenters.values(triangle_cell)
# logging.info("".join([" " for l in range(loglevel)])+triangle_epidermis_barycenters.shape
triangle_epidermis_directions = triangle_epidermis_barycenters - triangle_cell_barycenters
# logging.info("".join([" " for l in range(loglevel)])+triangle_epidermis_directions.shape
triangle_normal_vectors = topomesh.wisp_property('normal',degree=2).values(epidermis_triangles)
# reversed_normals = np.where(np.dot(triangle_normal_vectors,triangle_interface_directions) < 0)[0]
reversed_normals = np.where(np.einsum('ij,ij->i',triangle_normal_vectors,triangle_epidermis_directions) < 0)[0]
triangle_normal_vectors[reversed_normals] = -triangle_normal_vectors[reversed_normals]
# cell_normal_vectors = np.transpose([nd.mean(triangle_normal_vectors[:,0],triangle_cell,index=list(topomesh.wisps(3))),
# nd.mean(triangle_normal_vectors[:,1],triangle_cell,index=list(topomesh.wisps(3))),
# nd.mean(triangle_normal_vectors[:,2],triangle_cell,index=list(topomesh.wisps(3)))])
# cell_normal_vectors = cell_normal_vectors/np.linalg.norm(cell_normal_vectors,axis=1)[:,np.newaxis]
# cell_normal_vectors = array_dict(cell_normal_vectors,keys=list(topomesh.wisps(3)))
# triangle_target_normals = cell_normal_vectors.values(triangle_cell)
triangle_radius_direction = triangle_barycenters - triangle_cell_barycenters
triangle_target_normals = triangle_radius_direction/np.maximum(1e-5,np.linalg.norm(triangle_radius_direction,axis=1))[:,np.newaxis]
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(epidermis_triangles)
epidermis_vertices = np.unique(triangle_vertices)
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_point_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_vertices[:,0]) - np.tile(triangle_epidermis_barycenters,(3,1))
triangle_point_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_vertices[:,0]) - np.tile(triangle_barycenters,(3,1))
triangle_plane_products = np.einsum('ij,ij->i',triangle_point_vectors,np.tile(triangle_target_normals,(3,1)))
triangle_point_projectors = -triangle_plane_products[:,np.newaxis]*np.tile(triangle_target_normals,(3,1))
triangle_projectors_norm = np.linalg.norm(triangle_point_projectors,axis=1)
non_planar_points = np.where(triangle_projectors_norm > 1.)[0]
triangle_point_projectors[non_planar_points] = triangle_point_projectors[non_planar_points]/np.maximum(1e-5,triangle_projectors_norm)[non_planar_points,np.newaxis]
convex_points = np.where(triangle_plane_products>0.)[0]
triangle_point_projectors[convex_points] = np.array([0.,0.,0.])
convexity_unitary_force = triangle_point_projectors
epidermis_convexity_force = np.transpose([nd.sum(convexity_unitary_force[:,0],triangle_vertices[:,0],index=list(topomesh.wisps(0))),
nd.sum(convexity_unitary_force[:,1],triangle_vertices[:,0],index=list(topomesh.wisps(0))),
nd.sum(convexity_unitary_force[:,2],triangle_vertices[:,0],index=list(topomesh.wisps(0)))])
return epidermis_convexity_force
def _property_topomesh_epidermis_planarization_force(topomesh, verbose=False, debug=False, loglevel=0):
if not topomesh.has_wisp_property('epidermis',degree=2,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',degree=2)
if not topomesh.has_wisp_property('epidermis',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'epidermis',degree=3)
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('cells',degree=0,is_computed=True):
compute_topomesh_property(topomesh,'cells',degree=0)
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 not topomesh.has_wisp_property('length',degree=1,is_computed=True):
compute_topomesh_property(topomesh,'length',degree=1)
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)
epidermis_triangles = np.array(list(topomesh.wisps(2)))[topomesh.wisp_property('epidermis',2).values()]
epidermis_cells = np.array(list(topomesh.wisps(3)))[topomesh.wisp_property('epidermis',3).values()]
epidermis_vertices = np.array(list(topomesh.wisps(0)))[topomesh.wisp_property('epidermis',0).values()]
triangle_cells = np.concatenate(topomesh.wisp_property('cells',2).values(epidermis_triangles))
triangle_areas = topomesh.wisp_property('area',2).values(epidermis_triangles)
triangle_centers = topomesh.wisp_property('barycenter',2).values(epidermis_triangles)
triangle_weighted_normals = triangle_areas[:,np.newaxis]*topomesh.wisp_property('normal',2).values(epidermis_triangles)
cell_areas = nd.sum(triangle_areas,triangle_cells,index=epidermis_cells)
cell_epidermis_centers = np.transpose([nd.sum(triangle_areas*triangle_centers[:,k],triangle_cells,index=epidermis_cells) for k in range(3)])
cell_epidermis_centers = cell_epidermis_centers/cell_areas[:,np.newaxis]
cell_epidermis_centers = array_dict(cell_epidermis_centers,epidermis_cells)
cell_weighted_normals = np.transpose([nd.sum(triangle_weighted_normals[:,k],triangle_cells,index=epidermis_cells) for k in range(3)])
cell_normals = cell_weighted_normals/cell_areas[:,np.newaxis]
cell_normals = array_dict(cell_normals/np.linalg.norm(cell_normals,axis=1)[:,np.newaxis],epidermis_cells)
vertex_cells = np.concatenate([np.intersect1d(np.array(topomesh.wisp_property('cells',0)[v],int),epidermis_cells) for v in epidermis_vertices])
vertex_cell_vertex = np.concatenate([v*np.ones_like(np.intersect1d(topomesh.wisp_property('cells',0)[v],epidermis_cells),int) for v in epidermis_vertices])
vertex_cell_normals = cell_normals.values(vertex_cells)
#vertex_cell_normals = np.array([cell_normals[c] if cell_normals.has_key(c) else np.array([0,0,1]) for c in vertex_cells])
vertex_cell_centers = cell_epidermis_centers.values(vertex_cells)
#vertex_cell_centers = np.array([cell_epidermis_centers[c] if cell_epidermis_centers.has_key(c) else topomesh.wisp_property('barycenter',3)[c]])
vertex_cell_vectors = topomesh.wisp_property('barycenter',0).values(vertex_cell_vertex) - vertex_cell_centers
vertex_cell_projectors = -np.einsum('ij,ij->i',vertex_cell_vectors,vertex_cell_normals)[:,np.newaxis]*vertex_cell_normals
vertex_projectors = np.transpose([nd.sum(vertex_cell_projectors[:,k],vertex_cell_vertex,index=list(topomesh.wisps(0))) for k in range(3)])
vertex_projectors[np.isnan(vertex_projectors)] = 0.
return vertex_projectors
def _property_topomesh_cell_interface_planarization_force(topomesh, verbose=False, debug=False, loglevel=0):
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('barycenter',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'barycenter',degree=3)
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('neighbors',degree=3,is_computed=True):
compute_topomesh_property(topomesh,'neighbors',degree=3)
planarization_force = np.zeros_like(topomesh.wisp_property('barycenter',degree=0).values(),np.float32)
for c in topomesh.wisps(3):
logging.info("".join([" " for l in range(loglevel)])+"Cell "+str(c)+" interface planarization")
for n in topomesh.wisp_property("neighbors",degree=3)[c]:
normal_direction = topomesh.wisp_property("barycenter",degree=3)[n] - topomesh.wisp_property("barycenter",degree=3)[c]
# triangle_vertices = np.array([list(topomesh.borders(2,t,2)) for t in topomesh.borders(3,c) if t in topomesh.borders(3,n)])
interface_triangles = np.intersect1d(topomesh.wisp_property('borders',degree=3)[c],topomesh.wisp_property('borders',degree=3)[n])
triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(interface_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)
# interface_vertices = np.intersect1d(cell_vertices[c],cell_vertices[n])
interface_vertices = np.unique(triangle_vertices)
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_directions = triangle_edge_vectors/triangle_edge_lengths[...,np.newaxis]
normal_vectors = np.cross(triangle_edge_vectors[:,1],triangle_edge_vectors[:,2])
# reversed_normals = np.where(np.dot(normal_vectors,normal_direction) < 0)[0]
reversed_normals = np.where(np.core.umath_tests.inner1d(normal_vectors,normal_direction) < 0)[0]
normal_vectors[reversed_normals] = np.cross(triangle_edge_vectors[reversed_normals,2],triangle_edge_vectors[reversed_normals,1])
normal_vectors = normal_vectors/np.linalg.norm(normal_vectors,axis=1)[:,np.newaxis]
if(np.isnan(normal_vectors).any()):
logging.info("".join([" " for l in range(loglevel)])+"Interface "+str(c)+" -> "+str(n)+" planarization force")
normal_target = np.mean(normal_vectors,axis=0)
normal_target = normal_target/np.linalg.norm(normal_target)
if(np.isnan(normal_target).any()):
logging.info("".join([" " for l in range(loglevel)])+"Interface "+str(c)+" -> "+str(n)+" planarization force")
# normal_errors = 1.0 - np.dot(normal_vectors,normal_target)
barycenter = np.mean(topomesh.wisp_property('barycenter',degree=0).values(interface_vertices),axis=0)
point_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_vertices[:,0]) - barycenter
point_projectors = -np.dot(point_vectors,normal_target)[:,np.newaxis]*normal_target
non_planar_points = np.where(abs(np.dot(point_vectors,normal_target)) > 1.)[0]
point_projectors[non_planar_points] = point_projectors[non_planar_points]/abs(np.dot(point_vectors,normal_target))[non_planar_points,np.newaxis]
# planarization_unitary_force = normal_errors[:,np.newaxis]*point_projectors
planarization_unitary_force = point_projectors
if(np.isnan(planarization_unitary_force).any()):
logging.info("".join([" " for l in range(loglevel)])+"Interface "+str(c)+" -> "+str(n)+" planarization force")
interface_planarization_force = np.transpose([nd.sum(planarization_unitary_force[:,0],triangle_vertices[:,0],index=interface_vertices),
nd.sum(planarization_unitary_force[:,1],triangle_vertices[:,0],index=interface_vertices),
nd.sum(planarization_unitary_force[:,2],triangle_vertices[:,0],index=interface_vertices)])
planarization_force[interface_vertices] += interface_planarization_force
planarization_force[np.isnan(planarization_force)] = 0.
return planarization_force
[docs]def topomesh_triangle_split(input_topomesh, verbose=False, debug=False, loglevel=0):
"""Remesh a triangular PropertyTopomesh by spltting all faces into 4
Performs a simple remeshing of a triangular topomesh by inserting a middle
vertex on each edge, and splitting each triangle into four identical
triangle by linking the middle vertices.
>>> " • • "
>>> " / \ / \ "
>>> " / \ ---> •---• "
>>> " / \ / \ / \ "
>>> "•-------• •---•---•"
Parameters
----------
input_topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The initial structure on which to perform the remeshing
Returns
-------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The remeshed structure
Warnings
--------
The PropertyTopomesh must be a *triangular topomesh* (:meth:`cellcomplex.property_topomesh.analysis.is_triangular`)
"""
assert is_triangular(input_topomesh)
topomesh = deepcopy(input_topomesh)
triangle_vertex_positions = topomesh.wisp_property('barycenter',0).to_dict()
compute_topomesh_property(topomesh,'length',1)
compute_topomesh_property(topomesh,'vertices',1)
topomesh_edges = list(topomesh.wisps(1))
compute_topomesh_property(topomesh,'borders',2)
compute_topomesh_property(topomesh,'vertices',2)
topomesh_triangles = list(topomesh.wisps(2))
start_time = time()
logging.info("".join([" " for l in range(loglevel)])+"--> Splitting edges")
edge_splitted = {}
edge_middles = {}
for e in topomesh_edges:
edge_vertices = topomesh.wisp_property('vertices',1)[e]
mid = topomesh.add_wisp(0)
triangle_vertex_positions[mid] = topomesh.wisp_property('barycenter',0).values(edge_vertices).mean(axis=0)
topomesh.unlink(1,e,edge_vertices[1])
topomesh.link(1,e,mid)
eid = topomesh.add_wisp(1)
topomesh.link(1,eid,edge_vertices[1])
topomesh.link(1,eid,mid)
edge_splitted[e] = [e,eid]
edge_middles[e] = mid
end_time = time()
logging.info("".join([" " for l in range(loglevel)])+"--> Splitting edges ["+str(end_time-start_time)+" s]")
start_time = time()
logging.info("".join([" " for l in range(loglevel)])+"--> Splitting triangles")
for t in topomesh_triangles:
for e in topomesh.wisp_property('borders',2)[t]:
topomesh.unlink(2,t,e)
eid = topomesh.add_wisp(1)
for neighbor_edge in topomesh.wisp_property('borders',2)[t]:
if neighbor_edge != e:
topomesh.link(1,eid,edge_middles[neighbor_edge])
topomesh.link(2,t,eid)
fid = topomesh.add_wisp(2)
topomesh.link(2,fid,eid)
opposite_pid = list(set(topomesh.wisp_property('vertices',2)[t]).difference(topomesh.wisp_property('vertices',1)[e]))[0]
step_start_time = time()
for neighbor_edge in topomesh.wisp_property('borders',2)[t]:
if neighbor_edge != e:
for opposite_edge in edge_splitted[neighbor_edge]:
if opposite_pid in topomesh.borders(1,opposite_edge):
topomesh.link(2,fid,opposite_edge)
for c in topomesh.regions(2,t):
topomesh.link(3,c,fid)
end_time = time()
logging.info("".join([" " for l in range(loglevel)])+"--> Splitting triangles ["+str(end_time-start_time)+" s]")
topomesh.update_wisp_property('barycenter',degree=0,values=triangle_vertex_positions)
compute_topomesh_property(topomesh,'length',1)
compute_topomesh_property(topomesh,'vertices',1)
compute_topomesh_property(topomesh,'borders',2)
compute_topomesh_property(topomesh,'vertices',2)
return topomesh
[docs]def property_topomesh_edge_flip_optimization(topomesh,omega_energies=dict([('regularization',0.15),('neighborhood',0.65)]),simulated_annealing=True,display=False, image_edges=None, verbose=False, debug=False, loglevel=0,**kwargs):
"""Perform a topological optimization of a PropertyTopomesh using edge flips.
An energy functional is iteratively minimized by successively flipping edges
of the triangular topomesh. Edge flips that allow a decrease of the energy
while keeping the structure valid topologically and geometrically are
applied in ascending order of energy variation. To avoid local minima of
energy it is possible to apply energy increasing flips in a simulated
annealing procedure.
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to perform the optimization
omega_energies : dict
The weights (float) associated to each energy term:
* *regularization* energy optimizing triangle eccentricity
* *length* energy minimizing total edge length
* *neighborhood* energy optimizing vertex valency towards 6
* *image* energy optimizing adjacencies to observed ones
simulated_annealing : bool
Whether to allow or not energy increasing edge flips in the process
image_edges : :class:`numpy.ndarray`
The list of pairs of vertices used as reference if *image* energy is used
Returns
-------
n_flips : int
Number of edges flipped in the last iteration
Note
-------
The PropertyTopomesh passed as argument is updated.
Warnings
--------
The PropertyTopomesh must be a *triangular topomesh* (:meth:`cellcomplex.property_topomesh.analysis.is_triangular`)
"""
projected = kwargs.get("projected_map_display",False)
if image_edges is None:
img_edges = np.array([])
else:
img_edges = deepcopy(image_edges)
if simulated_annealing:
iterations = kwargs.get('iterations',20)
lambda_temperature = kwargs.get('lambda_temperature',0.85)
minimal_temperature = kwargs.get('minimal_temperature',0.1)
initial_temperature = minimal_temperature*np.power(1./lambda_temperature,iterations)
else:
iterations = kwargs.get('iterations',3)
initial_temperature = 1.
lambda_temperature = 0.9
minimal_temperature = np.power(lambda_temperature,iterations)
simulated_annealing_temperature = initial_temperature
while simulated_annealing_temperature > minimal_temperature:
simulated_annealing_temperature *= lambda_temperature
compute_topomesh_property(topomesh,'faces',0)
compute_topomesh_property(topomesh,'borders',1)
compute_topomesh_property(topomesh,'vertices',1)
compute_topomesh_property(topomesh,'faces',1)
compute_topomesh_property(topomesh,'length',1)
compute_topomesh_property(topomesh,'borders',2)
compute_topomesh_property(topomesh,'vertices',2)
compute_topomesh_triangle_properties(topomesh)
# logging.info("".join([" " for l in range(loglevel)])+"Area : ",topomesh.wisp_property('area',2).values().mean()," (",topomesh.wisp_property('area',2).values().std()," ) Eccentricity : ,",topomesh.wisp_property('eccentricity',2).values().mean()," [ T = ",simulated_annealing_temperature," ]"
flippable_edges = topomesh.wisp_property('faces',1).keys()[np.where(np.array(list(map(len,topomesh.wisp_property('faces',1).values()))) == 2)]
flippable_edge_vertices = topomesh.wisp_property('vertices',1).values(flippable_edges)
flippable_edge_triangle_vertices = topomesh.wisp_property('vertices',2).values(topomesh.wisp_property('faces',1).values(flippable_edges))
flippable_edge_flipped_vertices = np.array([list(set(list(np.unique(t))).difference(v)) for t,v in zip(flippable_edge_triangle_vertices,flippable_edge_vertices)])
wrong_edges = np.where(np.array(list(map(len,flippable_edge_flipped_vertices))) != 2)[0]
flippable_edges = np.delete(flippable_edges,wrong_edges,0)
flippable_edge_vertices = np.delete(flippable_edge_vertices,wrong_edges,0)
flippable_edge_triangle_vertices = np.delete(flippable_edge_triangle_vertices,wrong_edges,0)
flippable_edge_flipped_vertices = np.array([e for e in np.delete(flippable_edge_flipped_vertices,wrong_edges,0)])
flippable_edge_flipped_triangle_vertices = np.array([[np.concatenate([f,[v]]) for v in e] for (e,f) in zip(flippable_edge_vertices,flippable_edge_flipped_vertices)])
flippable_edge_triangle_areas = np.concatenate([triangle_geometric_features(flippable_edge_triangle_vertices[:,e],topomesh.wisp_property('barycenter',0),features=['area']) for e in [0,1]],axis=1)
flippable_edge_flipped_triangle_areas = np.concatenate([triangle_geometric_features(flippable_edge_flipped_triangle_vertices[:,e],topomesh.wisp_property('barycenter',0),features=['area']) for e in [0,1]],axis=1)
flippable_edge_flipped_triangle_areas[np.isnan(flippable_edge_flipped_triangle_areas)] = 100.
average_area = np.nanmean(topomesh.wisp_property('area',2).values())
flat = topomesh.wisp_property('barycenter',0).values()[:,2].std() < 1e-6
if flat:
wrong_edges = np.where(np.abs(flippable_edge_triangle_areas.sum(axis=1)-flippable_edge_flipped_triangle_areas.sum(axis=1)) > 1e-6)
else:
wrong_edges = np.where(np.abs(flippable_edge_triangle_areas.sum(axis=1)-flippable_edge_flipped_triangle_areas.sum(axis=1)) > average_area/10.)
flippable_edges = np.delete(flippable_edges,wrong_edges,0)
flippable_edge_vertices = np.delete(flippable_edge_vertices,wrong_edges,0)
flippable_edge_triangle_vertices = np.delete(flippable_edge_triangle_vertices,wrong_edges,0)
flippable_edge_flipped_vertices = np.delete(flippable_edge_flipped_vertices,wrong_edges,0)
flippable_edge_flipped_triangle_vertices = np.delete(flippable_edge_flipped_triangle_vertices,wrong_edges,0)
flippable_edge_triangle_areas = np.delete(flippable_edge_triangle_areas,wrong_edges,0)
flippable_edge_flipped_triangle_areas = np.delete(flippable_edge_flipped_triangle_areas,wrong_edges,0)
flippable_edge_energy_variation = np.zeros_like(flippable_edges,np.float)
if 'regularization' in omega_energies:
flippable_edge_area_error = np.power(flippable_edge_triangle_areas-average_area,2.0).sum(axis=1)
flippable_edge_area_flipped_error = np.power(np.maximum(flippable_edge_flipped_triangle_areas-average_area,0),2.0).sum(axis=1)
flippable_edge_area_energy_variation = array_dict(0.02*(flippable_edge_area_flipped_error-flippable_edge_area_error)/np.power(4.0*average_area,2.0),flippable_edges)
# flippable_edge_area_energy_variation = array_dict(flippable_edge_flipped_triangle_areas.sum(axis=1)-flippable_edge_triangle_areas.sum(axis=1),flippable_edges)
flippable_edge_triangle_eccentricities = np.concatenate([triangle_geometric_features(flippable_edge_triangle_vertices[:,e],topomesh.wisp_property('barycenter',0),features=['sinus_eccentricity']) for e in [0,1]],axis=1)
flippable_edge_flipped_triangle_eccentricities = np.concatenate([triangle_geometric_features(flippable_edge_flipped_triangle_vertices[:,e],topomesh.wisp_property('barycenter',0),features=['sinus_eccentricity']) for e in [0,1]],axis=1)
flippable_edge_eccentricity_energy_variation = array_dict(8.0*(flippable_edge_flipped_triangle_eccentricities.sum(axis=1)-flippable_edge_triangle_eccentricities.sum(axis=1)),flippable_edges)
flippable_edge_energy_variation += omega_energies['regularization']*flippable_edge_area_energy_variation.values(flippable_edges)
flippable_edge_energy_variation += omega_energies['regularization']*flippable_edge_eccentricity_energy_variation.values(flippable_edges)
if 'length' in omega_energies:
flippable_edge_lengths = np.linalg.norm(topomesh.wisp_property('barycenter',0).values(flippable_edge_vertices[:,1]) - topomesh.wisp_property('barycenter',0).values(flippable_edge_vertices[:,0]),axis=1)
flippable_edge_flipped_lengths = np.linalg.norm(topomesh.wisp_property('barycenter',0).values(flippable_edge_flipped_vertices[:,1]) - topomesh.wisp_property('barycenter',0).values(flippable_edge_flipped_vertices[:,0]),axis=1)
flippable_edge_length_energy_variation = array_dict(flippable_edge_flipped_lengths-flippable_edge_lengths,flippable_edges)
flippable_edge_energy_variation += omega_energies['length']*flippable_edge_length_energy_variation.values(flippable_edges)
if 'neighborhood' in omega_energies:
compute_topomesh_property(topomesh,'valency',0)
nested_mesh = kwargs.get("nested_mesh",False)
if nested_mesh:
compute_topomesh_property(topomesh,'cells',0)
compute_topomesh_property(topomesh,'epidermis',0)
target_neighborhood = array_dict((np.array(list(map(len,topomesh.wisp_property('cells',0).values()))) + topomesh.wisp_property('epidermis',0).values())*3,list(topomesh.wisps(0)))
else:
target_neighborhood = array_dict(np.ones_like(list(topomesh.wisps(0)))*6,list(topomesh.wisps(0)))
flippable_edge_neighborhood_error = np.power(np.abs(topomesh.wisp_property('valency',0).values(flippable_edge_vertices)-target_neighborhood.values(flippable_edge_vertices)),2.0).sum(axis=1)
flippable_edge_neighborhood_error += np.power(np.abs(topomesh.wisp_property('valency',0).values(flippable_edge_flipped_vertices)-target_neighborhood.values(flippable_edge_flipped_vertices)),2.0).sum(axis=1)
flippable_edge_neighborhood_flipped_error = np.power(np.abs(topomesh.wisp_property('valency',0).values(flippable_edge_vertices)-1-target_neighborhood.values(flippable_edge_vertices)),2.0).sum(axis=1)
flippable_edge_neighborhood_flipped_error += np.power(np.abs(topomesh.wisp_property('valency',0).values(flippable_edge_flipped_vertices)+1-target_neighborhood.values(flippable_edge_flipped_vertices)),2.0).sum(axis=1)
flippable_edge_neighborhood_energy_variation = array_dict(flippable_edge_neighborhood_flipped_error-flippable_edge_neighborhood_error,flippable_edges)
flippable_edge_energy_variation += omega_energies['neighborhood']*flippable_edge_neighborhood_energy_variation.values(flippable_edges)
if 'image' in omega_energies:
if img_edges is not None and len(img_edges)>0:
flippable_edge_energy_variation += omega_energies['image']*((vq(flippable_edge_vertices,img_edges)[1]==0).astype(float) - (vq(flippable_edge_flipped_vertices,img_edges)[1]==0 ).astype(float))
flippable_edge_energy_variation = array_dict(flippable_edge_energy_variation,flippable_edges)
if display:
pass
flippable_edge_sorted_energy_variation = array_dict(np.sort(flippable_edge_energy_variation.values()),flippable_edges[np.argsort(flippable_edge_energy_variation.values())])
start_time = time()
logging.info("".join([" " for l in range(loglevel)])+"--> Flipping mesh edges")
flipped_edges = 0
flipped_nonoptimal_edges = 0
unflippable_edges = set()
for eid in flippable_edge_sorted_energy_variation.keys():
flip_probability = np.exp(-flippable_edge_sorted_energy_variation[eid]/simulated_annealing_temperature)
if (simulated_annealing and (np.random.rand() < flip_probability)) or (1 < flip_probability):
if not eid in unflippable_edges:
# logging.info("".join([" " for l in range(loglevel)])+topomesh.wisp_property('vertices',1)[eid], flippable_edge_sorted_energy_variation[e]
# neighbor_edges = list(topomesh.region_neighbors(1,eid))
flipped_edges += 1
flipped_nonoptimal_edges += (flip_probability<1)
neighbor_edges = list(topomesh.border_neighbors(1,eid))
topomesh_flip_edge(topomesh,eid)
# for n_eid in neighbor_edges:
# unflippable_edges.append(n_eid)
unflippable_edges = unflippable_edges.union(set(neighbor_edges))
end_time = time()
logging.info("".join([" " for l in range(loglevel)])+" --> Flipped "+str(flipped_edges)+" edges ("+str(flipped_nonoptimal_edges)+" non-optimal) [ T = "+str(simulated_annealing_temperature)+"]")
logging.info("".join([" " for l in range(loglevel)])+"<-- Flipping mesh edges ["+str(end_time-start_time)+" s]")
return flipped_edges
[docs]def property_topomesh_edge_split_optimization(topomesh, maximal_length=None, split_adjacent_faces=True, iterations=1, verbose=False, debug=False, loglevel=0):
"""Perform a topological optimization of a PropertyTopomesh using edge splits.
The edges of the topomesh are iteratively splitted in order to become all
lesser than a target maximal length, while keeping the structure valid
topologically and geometrically. The edges are splitted in the descending
order of their lengths, until no edge is longer than *maximal_length* or
the maximal number of *iterations* is reached.
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to perform the optimization
maximal_length : float
The edge length above which the edges will be splitted
split_adjacent_faces : bool
Whether to divide the faces adjacent to the splitted edges
iterations : int
Number of optimization passes to perform
Returns
-------
n_splits : int
Number of edges splitted in the last iteration
Note
-------
The PropertyTopomesh passed as argument is updated.
Warnings
--------
If the PropertyTopomesh must is a *triangular topomesh* the *split_adjacent_faces* argument should be `True`
"""
compute_topomesh_property(topomesh,'vertices',1)
compute_topomesh_property(topomesh,'length',1)
if maximal_length is None:
target_length = np.percentile(topomesh.wisp_property('length',1).values(),70)
maximal_length = 4./3. * target_length
for iteration in range(iterations):
sorted_edge_length_edges = np.array(list(topomesh.wisps(1)))[np.argsort(-topomesh.wisp_property('length',1).values(list(topomesh.wisps(1))))]
sorted_edge_length_edges = sorted_edge_length_edges[topomesh.wisp_property('length',1).values(sorted_edge_length_edges) > maximal_length]
modified_edges = set()
n_splits = 0
for eid in sorted_edge_length_edges:
if not eid in modified_edges:
#logging.info("".join([" " for l in range(loglevel)])+" <-- Splitting edge ",eid," [",np.min(map(len,map(np.unique,[list(topomesh.borders(1,e)) for e in topomesh.wisps(1)]))),"]"
modified_edges = modified_edges.union(set(np.unique(np.concatenate([np.array(list(topomesh.borders(2,fid))) for fid in topomesh.regions(1,eid)]))))
topomesh_split_edge(topomesh,eid,split_adjacent_faces=split_adjacent_faces)
n_splits += 1
#logging.info("".join([" " for l in range(loglevel)])+" <-- Splitted edge ",eid," [",np.min(map(len,map(np.unique,[list(topomesh.borders(1,e)) for e in topomesh.wisps(1)]))),"]"
logging.info("".join([" " for l in range(loglevel)])+"--> Splitted "+str(n_splits)+" edges")
compute_topomesh_property(topomesh,'vertices',1)
compute_topomesh_property(topomesh,'length',1)
return n_splits
[docs]def property_topomesh_edge_collapse_optimization(topomesh, omega_energies=dict([('length',0.01),('error_quadrics',0.65)]), minimal_length=None, target_triangles=None, iterations=1, verbose=False, debug=False, loglevel=0):
"""Perform a topological optimization of a PropertyTopomesh using edge collapses.
The topomesh is decimated by successive edge collapse operations in order to
minimize an energy function, while reaching the double goal of eliminating
edges under a minimal length and reducing the total number of triangles. The
collapse operations are performed in the ascending order of energy variation
until either the target number of triangles or the maximal number of *iterations*
is reached.
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to perform the optimization
omega_energies : dict
The weights (float) associated to each energy term:
* *error_quadics* energy minimizing the decimation error
* *length* energy maximizing total edge length
minimal_length : float
The edge length under which the edges will be collapsed
target_triangles : int
The number of triangles to reach after decimation
iterations : int
Number of optimization passes to perform
Returns
-------
n_collapses : int
Number of edges collapsed in the last iteration
Note
-------
The PropertyTopomesh passed as argument is updated.
Warnings
--------
The PropertyTopomesh must be a *triangular topomesh* (:meth:`cellcomplex.property_topomesh.analysis.is_triangular`)
"""
compute_topomesh_property(topomesh,'vertices',1)
compute_topomesh_property(topomesh,'vertices',2)
compute_topomesh_property(topomesh,'length',1)
compute_topomesh_property(topomesh,'area',2)
if minimal_length is None:
target_length = np.percentile(topomesh.wisp_property('length',1).values(),70)
minimal_length = 3./4. * target_length
if target_triangles is None:
target_triangles = topomesh.nb_wisps(2)/4
if 'error_quadrics' in omega_energies.keys():
maximal_quadrics_error = np.power(minimal_length,2.)
compute_topomesh_property(topomesh,'faces',0)
vertices_positions = topomesh.wisp_property('barycenter',0).values(topomesh.wisp_property('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_vectors = normal_vectors/normal_norms[:,np.newaxis]
plane_d = -np.einsum('...ij,...ij->...i',normal_vectors,vertices_positions[:,0])
triangle_planes = np.concatenate([normal_vectors,plane_d[:,np.newaxis]],axis=1)
triangle_plane_quadrics = array_dict(np.einsum('...i,...j->...ij',triangle_planes,triangle_planes),list(topomesh.wisps(2)))
vertex_triangles = deepcopy(topomesh.wisp_property('faces',0))
# vertex_quadrics = triangle_plane_quadrics.values([t[0]+1 for t in vertex_triangles])
# homogeneous_positions = np.concatenate([topomesh.wisp_property('barycenter',0).values(),np.ones((topomesh.nb_wisps(0),1))],axis=1)
for iterations in range(iterations):
compute_topomesh_property(topomesh,'vertices',1)
compute_topomesh_property(topomesh,'length',1)
edge_energy_variation = np.zeros_like(list(topomesh.wisps(1)),np.float)
if 'length' in omega_energies.keys():
edge_energy_variation += omega_energies['length']*(topomesh.wisp_property('length',1).values(list(topomesh.wisps(1))))
if 'error_quadrics' in omega_energies.keys():
edge_vertices = topomesh.wisp_property('vertices',1).values()
edge_middles = topomesh.wisp_property('barycenter',0).values(edge_vertices).mean(axis=1)
edge_middles_homogeneous = np.concatenate([edge_middles,np.ones((topomesh.nb_wisps(1),1))],axis=1)
edge_vertex_faces = np.array([np.concatenate(vertex_triangles.values(v)) for v in edge_vertices])
edge_vertex_face_quadrics = np.array([triangle_plane_quadrics.values(t) for t in edge_vertex_faces])
edge_vertex_face_middles = np.array([ [e for t in e_v_t] for e,e_v_t in zip(edge_middles_homogeneous,edge_vertex_faces)])
edge_quadrics_errors = np.array([np.abs(np.einsum('...ij,...ij->...i',m,np.einsum('...ij,...j->...i',q,m))).sum() for q,m in zip(edge_vertex_face_quadrics,edge_vertex_face_middles)])
edge_quadrics_errors = array_dict(edge_quadrics_errors,list(topomesh.wisps(1)))
edge_energy_variation += omega_energies['error_quadrics']*(edge_quadrics_errors.values(list(topomesh.wisps(1))))
sorted_energy_variation_edges = np.array(list(topomesh.wisps(1)))[np.argsort(edge_energy_variation)]
if 'error_quadrics' in omega_energies.keys():
sorted_energy_variation_edges = sorted_energy_variation_edges[edge_quadrics_errors.values(sorted_energy_variation_edges) < maximal_quadrics_error]
sorted_energy_variation_edges = sorted_energy_variation_edges[topomesh.wisp_property('length',1).values(sorted_energy_variation_edges) < minimal_length]
# sorted_quadrics_errors_edges = sorted_quadrics_errors_edges[np.sort(edge_quadrics_errors.values()) < np.percentile(edge_quadrics_errors.values(),5)]
modified_edges = set()
n_collapses = 0
for eid in sorted_energy_variation_edges:
if not eid in modified_edges and topomesh.nb_wisps(2)>target_triangles:
# logging.info("".join([" " for l in range(loglevel)])+" <-- Collapsing edge ",eid," [",np.min(map(len,map(np.unique,[list(topomesh.borders(1,e)) for e in topomesh.wisps(1)]))),"]"
modified_edges = modified_edges.union(set(np.unique(list(topomesh.border_neighbors(1,eid))))).union({eid})
collapsed = topomesh_collapse_edge(topomesh,eid,manifold=False)
n_collapses += collapsed
# logging.info("".join([" " for l in range(loglevel)])+" <-- Collapsed edge ",eid," [",np.min(map(len,map(np.unique,[list(topomesh.borders(1,e)) for e in topomesh.wisps(1)]))),"]"
logging.info("".join([" " for l in range(loglevel)])+"--> Collapsed "+str(n_collapses)+" edges ["+str(np.min(list(map(len,list(map(np.unique,[list(topomesh.borders(1,e)) for e in topomesh.wisps(1)]))))))+"]")
return n_collapses
[docs]def property_topomesh_isotropic_remeshing(initial_topomesh, maximal_length=None, minimal_length=None, collapse=False, iterations=1, verbose=False, debug=False, loglevel=0):
"""Perform a global optimization of a triangular PropertyTopomesh
The topomesh is remeshed iteratively by applying only local topological
operations and vertex deformation. Based on target minimal and maximal
edge lengths, the mesh undergoes a split and collapse optimization,
then a pass of edge flip optimization regularizes its topology and a
Taubin smoothing its geometry. The process is repeated to perform a
global isotropic remeshing of the initial mesh.
Parameters
----------
input_topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The initial structure on which to perform the remeshing
maximal_length : float, *optional*:
The edge length above which the edges will be splitted
minimal_length : float, *optional*
The edge length under which the edges will be collapsed
collapse : bool
Whether to use or not the edge collapse optimization
iterations : int
Number of optimization passes to perform
Note
--------
If no edge length is set, the target edge length *l* is set to the 8th decile of
initial edge length and the heuristic setting *minimal_length* to 3/4 *l* and
*maximal_length* to 4/3 *l* is used.
Returns
-------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The remeshed structure
Warnings
--------
The PropertyTopomesh must be a *triangular topomesh* (:meth:`cellcomplex.property_topomesh.analysis.is_triangular`)
"""
topomesh = deepcopy(initial_topomesh)
n_flips = topomesh.nb_wisps(1)
n_splits = topomesh.nb_wisps(1)
n_collapses = topomesh.nb_wisps(1)
compute_topomesh_property(topomesh,'length',1)
target_length = np.percentile(topomesh.wisp_property('length',1).values(),80)
if maximal_length is None:
maximal_length = 4./3. * target_length
if minimal_length is None:
minimal_length = 3./4. * target_length
iteration = 0
while (n_flips+n_splits+n_collapses > topomesh.nb_wisps(1)/100.) and (iteration<iterations):
if collapse:
n_collapses = property_topomesh_edge_collapse_optimization(topomesh, minimal_length=minimal_length, iterations=1)
else:
n_collapses = 0
n_splits = property_topomesh_edge_split_optimization(topomesh, maximal_length=maximal_length, iterations=1)
n_flips = property_topomesh_edge_flip_optimization(topomesh,omega_energies=dict([('neighborhood',0.65)]),simulated_annealing=False,iterations=5)
property_topomesh_vertices_deformation(topomesh,omega_forces=dict([('taubin_smoothing',0.33)]),iterations=5)
iteration += 1
return topomesh
