# -*- 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 Topological Operations
=============================
.. autosummary::
:toctree: generated/
topomesh_flip_edge -- Flip an edge element in a triangular PropertyTopomesh
topomesh_split_edge -- Split an edge element in a PropertyTopomesh
topomesh_collapse_edge -- Collapse an edge element in a triangular PropertyTopomesh
topomesh_split_vertex -- Split a vertex element in a triangular PropertyTopomesh
"""
import numpy as np
from scipy import ndimage as nd
# from scipy.cluster.vq import vq
from cellcomplex.utils.array_dict import array_dict
from cellcomplex.property_topomesh import PropertyTopomesh
from cellcomplex.property_topomesh.analysis import compute_topomesh_property, is_triangular
# from cellcomplex.property_topomesh.utils.tissue_analysis_tools import cell_vertex_extraction
from cellcomplex.property_topomesh.utils.array_tools import array_unique
# from vplants.container.topomesh_algo import is_collapse_topo_allowed, collapse_edge
from cellcomplex.property_topomesh.utils.geometry_tools import triangle_geometric_features
from time import time
from copy import deepcopy
import logging
[docs]def topomesh_remove_vertex(topomesh,pid,kept_fid=None, triangulate=True, verbose=False, debug=False, loglevel=0):
try:
vertex_fids = list(topomesh.regions(0,pid,2))
vertex_fid_cells = np.array([list(topomesh.regions(2,fid)) for fid in vertex_fids])
if len(vertex_fids) > 0:
assert array_unique(vertex_fid_cells).shape[0] == 1
if kept_fid is None:
kept_fid = vertex_fids[0]
fid_to_keep = kept_fid
fids_to_delete = list(set(vertex_fids).difference({fid_to_keep}))
vertex_eids = list(topomesh.regions(0,pid))
vertex_fid_eids = np.unique([list(topomesh.borders(2,fid)) for fid in vertex_fids])
eids_to_keep = list(set(vertex_fid_eids).difference(set(vertex_eids)))
for eid in vertex_eids:
topomesh.remove_wisp(1,eid)
for fid in fids_to_delete:
topomesh.remove_wisp(2,fid)
for eid in eids_to_keep:
if not eid in topomesh.borders(2,fid_to_keep):
topomesh.link(2,fid_to_keep,eid)
topomesh.remove_wisp(0,pid)
del topomesh.wisp_property("barycenter",0)[pid]
return True
except AssertionError:
logging.info("".join([" " for l in range(loglevel)])+"--> Impossible to remove vertex : wrong configuration ( ",array_unique(vertex_fid_cells).shape[0]," cell memberships)")
return False
[docs]def topomesh_collapse_edge(topomesh, eid, kept_pid=None, manifold=True, verbose=False, debug=False, loglevel=0):
"""Perform an edge collapse operation on a triangular PropertyTopomesh
The edge collapse operation merges the two vertices linked by a given mesh
edge, suppressed the incident faces and links the remaining elements to
preserve the topological consistency of the complex. The operation is
performed only if the resulting mesh is valid.
>>> •-----•-----• •-----•-----•
>>> |\ / \ /| | \ | / |
>>> | \ / \ / | | \ | / |
>>> •--•-----•--• --> •-----•-----•
>>> | / \ / \ | | / | \ |
>>> |/ \ / \| | / | \ |
>>> •-----•-----• •-----•-----•
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to perform the operation
eid : unsigned int
The ID of the edge element to collapse
kept_pid : unsigned int, *optional*
The ID of the edge vertex that should be kept after collapse
manifold : bool
Whether to consider the topomesh as a 2-manifold complex
Returns
-------
None
Note
-------
The PropertyTopomesh passed as argument is updated.
Warnings
--------
The PropertyTopomesh must be a *triangular topomesh* (:meth:`cellcomplex.property_topomesh.analysis.is_triangular`)
"""
initial_topomesh = deepcopy(topomesh)
try:
if manifold:
assert len(list(topomesh.regions(1,eid))) == 2
pid_to_keep, pid_to_delete = topomesh.borders(1,eid)
logging.info("".join([" " for l in range(loglevel)])+"--> Trying to collapse edge "+str(eid)+" : "+str(pid_to_keep)+" ; "+str(pid_to_delete))
if kept_pid is not None and pid_to_keep != kept_pid:
pid_to_delete, pid_to_keep = topomesh.borders(1,eid)
edge_fids = list(topomesh.regions(1,eid))
eids_to_keep = list(topomesh.regions(0,pid_to_keep))
eids_to_delete = list(topomesh.regions(0,pid_to_delete))
# logging.info("".join([" " for l in range(loglevel)])+" --> Edge fids : ",edge_fids)
assert np.all([len(list(set(topomesh.borders(2,fid)).intersection(set(eids_to_keep)))) == 2 for fid in edge_fids])
assert np.all([len(list(set(topomesh.borders(2,fid)).intersection(set(eids_to_delete)))) == 2 for fid in edge_fids])
if manifold:
assert np.all(np.concatenate([[len(list(topomesh.regions(1,e))) == 2 for e in topomesh.borders(2,fid)] for fid in edge_fids]))
pids_to_link = np.unique([list(set(topomesh.borders(2,fid,2)).difference({pid_to_keep,pid_to_delete}))[0] for fid in edge_fids])
# logging.info("".join([" " for l in range(loglevel)])+" --> Face pids : ",pids_to_link)
assert len(pids_to_link) == len(edge_fids)
# assert np.all(np.array(map(len,[list(topomesh.regions(0,pid)) for pid in pids_to_link])) > 3)
assert np.all([np.all([len(list(topomesh.regions(1,e))) == 2 for e in topomesh.borders(2,fid) if e != eid]) for fid in edge_fids])
for fid in edge_fids:
# logging.info("".join([" " for l in range(loglevel)])+" --> Face ",fid," : ",list(topomesh.borders(2,fid))," (",eids_to_keep,")"
eid_to_keep = list(set(topomesh.borders(2,fid)).intersection(set(eids_to_keep)).difference({eid}))[0]
eid_to_delete = list(set(topomesh.borders(2,fid)).intersection(set(eids_to_delete)).difference({eid}))[0]
# logging.info("".join([" " for l in range(loglevel)])+" --> Kept eid : ",eid_to_keep,list(topomesh.regions(1,eid_to_keep))
# logging.info("".join([" " for l in range(loglevel)])+" --> Deleted eid : ",eid_to_delete,list(topomesh.regions(1,eid_to_delete))
topomesh.unlink(2,fid,eid)
topomesh.unlink(2,fid,eid_to_keep)
topomesh.unlink(2,fid,eid_to_delete)
fid_to_link = list(topomesh.regions(1,eid_to_delete))[0]
#for fid_to_link in topomesh.regions(1,eid_to_delete):
topomesh.unlink(2,fid_to_link,eid_to_delete)
topomesh.link(2,fid_to_link,eid_to_keep)
topomesh.remove_wisp(1,eid_to_delete)
topomesh.remove_wisp(2,fid)
for eid_to_link in set(topomesh.regions(0,pid_to_delete)).difference({eid}):
topomesh.unlink(1,eid_to_link,pid_to_delete)
neighbor_pid = list(topomesh.borders(1,eid_to_link))[0]
if neighbor_pid not in topomesh.region_neighbors(0,pid_to_keep):
# topomesh.unlink(1,eid_to_link,pid_to_delete)
topomesh.link(1,eid_to_link,pid_to_keep)
else:
for fid_to_link in topomesh.regions(1,eid_to_link):
if pid_to_keep in topomesh.borders(2,fid_to_link,2):
topomesh.remove_wisp(2,fid_to_link)
else:
eid_to_fuse = tuple(set(topomesh.regions(0,pid_to_keep)).intersection(set(topomesh.regions(0,neighbor_pid))))[0]
topomesh.unlink(2,fid_to_link,eid_to_link)
topomesh.link(2,fid_to_link,eid_to_fuse)
topomesh.remove_wisp(1,eid_to_link)
if kept_pid is not None:
topomesh.wisp_property("barycenter",0)[pid_to_keep] = topomesh.wisp_property("barycenter",0)[kept_pid]
else:
topomesh.wisp_property("barycenter",0)[pid_to_keep] = (topomesh.wisp_property("barycenter",0).values([pid_to_keep,pid_to_delete]).sum(axis=0))/2.
topomesh.remove_wisp(0,pid_to_delete)
del topomesh.wisp_property("barycenter",0)[pid_to_delete]
topomesh.remove_wisp(1,eid)
# if np.max([topomesh.nb_regions(1,e) for e in topomesh.wisps(1)])>2:
# logging.info("".join([" " for l in range(loglevel)])+" --> Error while collapsing! Non-manifold edges :",np.array(list(topomesh.wisps(1)))[np.array([topomesh.nb_regions(1,e) for e in topomesh.wisps(1)])>2]
assert np.max([topomesh.nb_regions(1,e) for e in topomesh.wisps(1)])==2
logging.info("".join([" " for l in range(loglevel)])+"<-- Collapsed edge"+str(eid)+" : "+str(pid_to_keep)+" ("+str(topomesh.nb_wisps(2))+" Faces )")
#raw_input()
return True
# edge_vertices = np.sort(np.array([list(topomesh.borders(1,e)) for e in topomesh.wisps(1) if topomesh.nb_borders(1,e) == 2]))
# edge_vertex_id = edge_vertices[:,0]*10000 + edge_vertices[:,1]
# if edge_vertices.shape[0] != array_unique(edge_vertices).shape[0]:
# logging.info("".join([" " for l in range(loglevel)])+eid," collapse error : (",pid_to_keep,pid_to_delete,@")",np.array(list(topomesh.wisps(1)))[nd.sum(np.ones_like(edge_vertex_id),edge_vertex_id,index=edge_vertex_id)>1]
# raw_input()
except AssertionError:
logging.info("".join([" " for l in range(loglevel)])+"<-- Impossible to collapse edge : wrong configuration ( "+str(len(list(initial_topomesh.regions(1,eid))))+" regions)")
assert np.max([initial_topomesh.nb_regions(1,e) for e in initial_topomesh.wisps(1)])==2
topomesh._borders = deepcopy(initial_topomesh._borders)
topomesh._regions = deepcopy(initial_topomesh._regions)
topomesh.update_wisp_property('barycenter',0,initial_topomesh.wisp_property('barycenter',0))
assert np.max([topomesh.nb_regions(1,e) for e in topomesh.wisps(1)])==2
logging.info("".join([" " for l in range(loglevel)])+"<-- Failed to collapse edge"+str(eid)+" : "+str(pid_to_keep)+" ("+str(topomesh.nb_wisps(2))+" Faces )")
# raw_input()
return False
# def topomesh_collapse_edge(topomesh,eid,kept_pid=None):
# try:
# assert len(list(topomesh.regions(1,eid))) == 2
# pid_to_keep, pid_to_delete = topomesh.borders(1,eid)
# logging.info("".join([" " for l in range(loglevel)])+"Collapsing edge eid : ",pid_to_keep, pid_to_delete
# if kept_pid is not None and pid_to_keep != kept_pid:
# pid_to_delete, pid_to_keep = topomesh.borders(1,eid)
# edge_fids = list(topomesh.regions(1,eid))
# for fid in edge_fids:
# topomesh.remove_wisp(2,fid)
# for pid in [pid_to_keep, pid_to_delete]:
# topomesh.unlink(1,eid,pid)
# neighbor_edges = list(topomesh.regions(0,pid_to_delete))
# for neighbor_eid in neighbor_edges:
# topomesh.unlink(1,neighbor_eid,pid_to_delete)
# neighbor_pid = list(topomesh.borders(1,neighbor_eid))[0]
# logging.info("".join([" " for l in range(loglevel)])+neighbor_eid," : ", neighbor_pid, list(topomesh.region_neighbors(0,pid_to_keep))
# if neighbor_pid not in topomesh.region_neighbors(0,pid_to_keep):
# topomesh.link(1,neighbor_eid,pid_to_keep)
# logging.info("".join([" " for l in range(loglevel)])+neighbor_eid, list(topomesh.borders(1,neighbor_eid)),"[",list(topomesh.regions(1,neighbor_eid)),"]"
# else:
# eid_to_fuse = tuple(set(topomesh.regions(0,pid_to_keep)).intersection(set(topomesh.regions(0,neighbor_pid))))[0]
# logging.info("".join([" " for l in range(loglevel)])+neighbor_eid," -> ",tuple(set(topomesh.regions(0,pid_to_keep)).intersection(set(topomesh.regions(0,neighbor_pid))))
# topomesh.unlink(1,neighbor_eid,neighbor_pid)
# for fid in topomesh.regions(1,neighbor_eid):
# topomesh.unlink(2,fid,neighbor_eid)
# topomesh.link(2,fid,eid_to_fuse)
# topomesh.remove_wisp(1,neighbor_eid)
# logging.info("".join([" " for l in range(loglevel)])+eid_to_fuse, list(topomesh.borders(1,eid_to_fuse)),"[",list(topomesh.regions(1,eid_to_fuse)),"]"
# topomesh.wisp_property("barycenter",0)[pid_to_keep] = (topomesh.wisp_property("barycenter",0).values([pid_to_keep,pid_to_delete]).sum(axis=0))/2.
# topomesh.remove_wisp(0,pid_to_delete)
# del topomesh.wisp_property("barycenter",0)[pid_to_delete]
# topomesh.remove_wisp(1,eid)
# edge_borders = np.array(map(len,map(np.unique,[list(topomesh.borders(1,e)) for e in topomesh.wisps(1)])))
# if np.min(edge_borders) == 1:
# logging.info("".join([" " for l in range(loglevel)])+eid," collapse error - borders : (",pid_to_keep,pid_to_delete,")",np.array(list(topomesh.wisps(1)))[edge_borders==1],neighbor_edges
# raw_input()
# edge_regions = np.array(map(len,map(np.unique,[list(topomesh.regions(1,e)) for e in topomesh.wisps(1)])))
# if np.max(edge_regions) > 2:
# logging.info("".join([" " for l in range(loglevel)])+eid," collapse error - regions : (",pid_to_keep,pid_to_delete,")",np.array(list(topomesh.wisps(1)))[edge_regions>2],neighbor_edges
# raw_input()
# except AssertionError:
# logging.info("".join([" " for l in range(loglevel)])+"Impossible to collapse edge : wrong configuration ( ",len(list(topomesh.regions(1,eid)))," regions)"
# return False
[docs]def topomesh_flip_edge(topomesh, eid, verbose=False, debug=False, loglevel=0):
"""Perform an edge flip operation on a triangular PropertyTopomesh
The edge flip operation changes the vertices of a mesh edge to the two
opposite vertices of its adjacent triangles. Thel local topological
configuration has to be valid for the operation to be performed.
>>> •---• •---•
>>> |\ | | /|
>>> | \ | --> | / |
>>> | \| |/ |
>>> •---• •---•
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to perform the operation
eid : unsigned int
The ID of the edge element to flip
Returns
-------
None
Note
-------
The PropertyTopomesh passed as argument is updated.
Warnings
--------
The PropertyTopomesh must be a *triangular topomesh* (:meth:`cellcomplex.property_topomesh.analysis.is_triangular`)
"""
try:
assert len(list(topomesh.regions(1,eid))) == 2
edge_triangles = list(topomesh.regions(1,eid))
assert set(list(topomesh.regions(2,edge_triangles[0]))) == set(list(topomesh.regions(2,edge_triangles[1])))
edge_vertices = np.array(list(topomesh.borders(1,eid)))
edge_triangle_edges = list(topomesh.region_neighbors(1,eid))
edge_triangle_edge_vertices = [list(topomesh.borders(1,e)) for e in edge_triangle_edges]
new_triangle_edges = {}
for i_pid,pid in enumerate(edge_vertices):
new_triangle_edges[i_pid] = []
for e in edge_triangle_edges:
if pid in topomesh.borders(1,e):
new_triangle_edges[i_pid].append(e)
edge_triangle_vertices = np.unique(edge_triangle_edge_vertices)
flipped_edge_vertices = [pid for pid in edge_triangle_vertices if not pid in topomesh.borders(1,eid)]
assert len(flipped_edge_vertices) == 2
for pid in topomesh.borders(1,eid):
topomesh.unlink(1,eid,pid)
for pid in flipped_edge_vertices:
topomesh.link(1,eid,pid)
for i_pid,fid in enumerate(topomesh.regions(1,eid)):
for e in topomesh.borders(2,fid):
if e != eid:
topomesh.unlink(2,fid,e)
for e in new_triangle_edges[i_pid]:
topomesh.link(2,fid,e)
return True
except AssertionError:
logging.info("".join([" " for l in range(loglevel)])+"--> Impossible to flip edge : wrong configuration ( "+str(len(list(topomesh.regions(1,eid))))+" Faces)")
return False
[docs]def topomesh_split_edge(topomesh, eid, split_adjacent_faces=True, verbose=False, debug=False, loglevel=0):
"""Perform an edge split operation on a PropertyTopomesh
The edge split operation creates a new vertex in the middle of the
splitted edge, and creates two new edges to link the middle vertex to
the opposite vertices of the adjacent triangles. If the adjacent faces
are not splitted, the mesh is no longer triangular. On the other hand,
if the mesh is polygonal the splitting of the faces is not unique.
>>> •---• •---•
>>> |\ | |\ /|
>>> | \ | --> | • |
>>> | \| |/ \|
>>> •---• •---•
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to perform the operation
eid : unsigned int
The ID of the edge element to split
split_adjacent_faces : bool
Whether to split the adjacent face elements
Returns
-------
None
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`
"""
pid_to_keep, pid_to_unlink = topomesh.borders(1,eid)
edge_fids = list(topomesh.regions(1, eid))
if split_adjacent_faces:
eids_to_unlink = np.array([list(set(list(topomesh.borders(2,fid))).intersection(set(list(topomesh.regions(0,pid_to_unlink)))).difference({eid}))[0] for fid in edge_fids])
pids_split = np.array([list(set(list(topomesh.borders(2,fid,2))).difference({pid_to_keep, pid_to_unlink}))[0] for fid in edge_fids])
pid_to_add = topomesh.add_wisp(0, np.max(list(topomesh.wisps(0)))+1)
topomesh.unlink(1,eid,pid_to_unlink)
topomesh.link(1,eid,pid_to_add)
logging.info("".join([" " for l in range(loglevel)])+" Edge "+str(eid)+" : "+str(pid_to_keep)+", "+str(pid_to_add))
topomesh.wisp_property("barycenter",0)[pid_to_add] = (topomesh.wisp_property("barycenter",0).values([pid_to_keep,pid_to_unlink]).sum(axis=0))/2.
eid_to_add = topomesh.add_wisp(1,np.max(list(topomesh.wisps(1)))+1)
topomesh.link(1,eid_to_add,pid_to_add)
topomesh.link(1,eid_to_add,pid_to_unlink)
logging.info("".join([" " for l in range(loglevel)])+" Split "+str(eid_to_add)+" : "+str(pid_to_add)+", "+str(pid_to_unlink))
if split_adjacent_faces:
for fid, eid_to_unlink, pid_split in zip(edge_fids,eids_to_unlink,pids_split):
eid_split = topomesh.add_wisp(1)
topomesh.link(1,eid_split,pid_to_add)
topomesh.link(1,eid_split,pid_split)
# logging.info("".join([" " for l in range(loglevel)])+"Added ",eid_split," : ",pid_to_add,pid_split
topomesh.unlink(2,fid,eid_to_unlink)
topomesh.link(2,fid,eid_split)
fid_to_add = topomesh.add_wisp(2)
topomesh.link(2,fid_to_add,eid_to_add)
topomesh.link(2,fid_to_add,eid_to_unlink)
topomesh.link(2,fid_to_add,eid_split)
for cid in topomesh.regions(2,fid):
topomesh.link(3,cid,fid_to_add)
else:
for fid in edge_fids:
topomesh.link(2,fid,eid_to_add)
# edge_borders = np.array(map(len,map(np.unique,[list(topomesh.borders(1,e)) for e in topomesh.wisps(1)])))
# if np.min(edge_borders) == 1:
# edge_vertices = np.sort(np.array([list(topomesh.borders(1,e)) for e in topomesh.wisps(1) if topomesh.nb_borders(1,e) == 2]))
# edge_vertex_id = edge_vertices[:,0]*10000 + edge_vertices[:,1]
# if edge_vertices.shape[0] != array_unique(edge_vertices).shape[0]:
# logging.info("".join([" " for l in range(loglevel)])+eid," split error : (",pid_to_keep,pid_to_unlink,")",np.array(list(topomesh.wisps(1)))[nd.sum(np.ones_like(edge_vertex_id),edge_vertex_id,index=edge_vertex_id)>1]
# raw_input()
return True
[docs]def topomesh_split_triangle(topomesh,fid, verbose=False, debug=False, loglevel=0):
assert(is_triangular(topomesh))
eid_to_keep, eid_to_unlink_1, eid_to_unlink_2 = topomesh.borders(2,fid)
face_pids = list(topomesh.borders(2,fid,2))
center_pid = topomesh.add_wisp(0)
pid_eids = {}
for pid in face_pids:
eid = topomesh.add_wisp(1)
topomesh.link(1,eid,pid)
topomesh.link(1,eid,center_pid)
pid_eids[pid] = eid
added_fids = []
for eid_to_unlink in [eid_to_unlink_1,eid_to_unlink_2]:
topomesh.unlink(2,fid,eid_to_unlink)
edge_fid = topomesh.add_wisp(2)
topomesh.link(2,edge_fid,eid_to_unlink)
for pid in topomesh.borders(1,eid_to_unlink):
topomesh.link(2,edge_fid,pid_eids[pid])
added_fids += [edge_fid]
for pid in topomesh.borders(1,eid_to_keep):
topomesh.link(2,fid,pid_eids[pid])
for fid_to_add in added_fids:
for cid in topomesh.regions(2,fid):
topomesh.link(3,cid,fid_to_add)
for property_name in topomesh.wisp_property_names(0):
# for property_name in ['barycenter']:
face_properties = topomesh.wisp_property(property_name,0).values(face_pids)
if ((face_properties.ndim > 1) and (len(np.unique(list(map(len,face_properties))))==1)) or (face_properties.dtype != np.object):
topomesh.wisp_property(property_name,0)[center_pid] = topomesh.wisp_property(property_name,0).values(face_pids).mean(axis=0)
return True
[docs]def topomesh_split_vertex(topomesh, pid, new_pid=None, splitting_axis=None, splitting_distance=None):
"""Perform an vertex split operation on a triangular PropertyTopomesh
The vertex split operation replaces a mesh vertex by a pair of vertices
linked by an edge, and creates the necessary higher degree elements to
ensure to topological consistency of the complex. The split occurs in
a given direction, with a given distance, that can both be estimated
automatically given the local configuration of the mesh.
>>> •-----•-----• •-----•-----•
>>> | \ | / | |\ / \ /|
>>> | \ | / | | \ / \ / |
>>> •-----•-----• --> •--•-----•--•
>>> | / | \ | | / \ / \ |
>>> | / | \ | |/ \ / \|
>>> •-----•-----• •-----•-----•
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The structure on which to perform the operation
pid : unsigned int
The ID of the vertex element to split
new_pid : unsigned int, *optional*
The ID of the new vertex inserted by the split
splitting_axis : :class:`numpy.ndarray`, *optional*
The 3D direction in which to perform the split
splitting_distance : :float, *optional*
The length of the edge element resulting from the split
Returns
-------
None
Note
-------
The PropertyTopomesh passed as argument is updated.
Warnings
--------
The PropertyTopomesh must be a *triangular topomesh* (:meth:`cellcomplex.property_topomesh.analysis.is_triangular`)
"""
assert is_triangular(topomesh)
neighbor_pids = np.array(list(topomesh.region_neighbors(0,pid)))
if splitting_axis is None:
neighbor_points = topomesh.wisp_property('barycenter',0).values(neighbor_pids)
eval,evec = np.linalg.eigh(np.cov(neighbor_points,rowvar=False))
splitting_axis = evec[:,np.argmax(np.abs(eval))]
splitting_axis = splitting_axis/np.linalg.norm(splitting_axis)
if splitting_distance is None:
neighbor_points = topomesh.wisp_property('barycenter',0).values(neighbor_pids)
neighbor_distance = np.linalg.norm(neighbor_points - topomesh.wisp_property('barycenter',0)[pid],axis=1)
splitting_distance = neighbor_distance.mean()/2.
point_normal = None
if topomesh.has_wisp_property('normal',0):
point_normal = topomesh.wisp_property('normal',0)[pid]
else:
neighbor_points = topomesh.wisp_property('barycenter',0).values(neighbor_pids)
eval, evec = np.linalg.eigh(np.cov(neighbor_points, rowvar=False))
point_normal = evec[:,np.argmin(np.abs(eval))]
splitting_orthogonal_axis = np.cross(splitting_axis,point_normal)
neighbor_vectors = topomesh.wisp_property('barycenter',0).values(neighbor_pids) - topomesh.wisp_property('barycenter',0)[pid]
neighbor_splitting_dot_products = np.einsum("...ij,...ij->...i", neighbor_vectors, [splitting_axis])
neighbor_splitting_side = dict(zip(neighbor_pids,np.sign(neighbor_splitting_dot_products)))
neighbor_splitting_ortho_side = np.sign(np.einsum("...ij,...ij->...i", neighbor_vectors, [splitting_orthogonal_axis]))
if (neighbor_splitting_ortho_side==-1).sum()>0:
left_neighbor = neighbor_pids[neighbor_splitting_ortho_side==-1][np.argmin(np.abs(neighbor_splitting_dot_products[neighbor_splitting_ortho_side==-1]))]
neighbor_splitting_side[left_neighbor]=0
if (neighbor_splitting_ortho_side==1).sum()>0:
right_neighbor = neighbor_pids[neighbor_splitting_ortho_side==1][np.argmin(np.abs(neighbor_splitting_dot_products[neighbor_splitting_ortho_side==1]))]
neighbor_splitting_side[right_neighbor]=0
pid_to_add = topomesh.add_wisp(0,new_pid)
topomesh.wisp_property('barycenter',0)[pid_to_add] = topomesh.wisp_property('barycenter',0)[pid] - splitting_distance*splitting_axis/2.
topomesh.wisp_property('barycenter',0)[pid] = topomesh.wisp_property('barycenter',0)[pid] + splitting_distance*splitting_axis/2.
neighbor_splitting_side[pid]=0
neighbor_splitting_side[pid_to_add]=0
eid_to_add = topomesh.add_wisp(1)
topomesh.link(1,eid_to_add,pid)
topomesh.link(1,eid_to_add,pid_to_add)
cids = topomesh.regions(0,pid,3)
for nid in neighbor_pids:
neid = list(set(list(topomesh.regions(0, pid))).intersection(set(list(topomesh.regions(0, nid)))))[0]
if neighbor_splitting_side[nid] == 0:
print(nid,"Side vertex [",neid,"]")
neid_to_add = topomesh.add_wisp(1)
topomesh.link(1, neid_to_add, nid)
topomesh.link(1, neid_to_add, pid_to_add)
for nfid in topomesh.regions(1, neid):
face_side = np.sum([neighbor_splitting_side[p] for p in topomesh.borders(2, nfid, 2)])
if face_side<0:
topomesh.unlink(2,nfid,neid)
topomesh.link(2,nfid,neid_to_add)
fid_to_add = topomesh.add_wisp(2)
topomesh.link(2, fid_to_add, eid_to_add)
topomesh.link(2, fid_to_add, neid)
topomesh.link(2, fid_to_add, neid_to_add)
for cid in cids:
topomesh.link(3, cid, fid_to_add)
for nid in neighbor_pids:
neid = list(set(list(topomesh.regions(0,pid))).intersection(set(list(topomesh.regions(0,nid)))))[0]
if neighbor_splitting_side[nid] == 1:
print(nid,"Positive vertex [",neid,"]")
pass
elif neighbor_splitting_side[nid] == -1:
print(nid,"Negative vertex [",neid,"]")
topomesh.unlink(1,neid,pid)
topomesh.link(1,neid,pid_to_add)
return pid_to_add
[docs]def topomesh_remove_interface_vertex(topomesh, pid, verbose=False, debug=False, loglevel=0):
try:
vertex_fids = list(topomesh.regions(0,pid,2))
if len(vertex_fids) > 0:
fid_to_keep = np.min(vertex_fids)
vertex_fid_cells = [list(topomesh.regions(2,fid)) for fid in vertex_fids]
vertex_cells = np.unique(vertex_fid_cells)
assert np.all([set(list(cids)) == set(list(vertex_cells)) for cids in vertex_fid_cells])
face_eids = np.concatenate([[eid for eid in topomesh.borders(2,fid) if not pid in topomesh.borders(1,eid)] for fid in vertex_fids])
oriented_face_eids = [face_eids[0]]
oriented_face_eid_orientations = [1]
candidate_eids = face_eids
while len(oriented_face_eids) < len(face_eids) and (len(candidate_eids) > 0):
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 = list(set(list(topomesh.regions(0,end_pid))).intersection(set(list(face_eids))).difference({current_eid}))
if len(candidate_eids)>0:
next_eid = candidate_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]
assert len(oriented_face_eids) == len(face_eids)
eids_to_remove = []
for fid in np.sort(vertex_fids):
for eid in topomesh.borders(2,fid):
topomesh.unlink(2,fid,eid)
if pid in topomesh.borders(1,eid):
eids_to_remove += [eid]
else:
topomesh.link(2,fid_to_keep,eid)
else:
eids_to_remove = list(topomesh.regions(0,pid))
# logging.info("".join([" " for l in range(loglevel)])+np.unique(eids_to_remove)
for eid in np.unique(eids_to_remove):
for pid_to_unlink in topomesh.borders(1,eid):
topomesh.unlink(1,eid,pid_to_unlink)
topomesh.remove_wisp(1,eid)
for fid in vertex_fids:
if fid != fid_to_keep:
topomesh.remove_wisp(2,fid)
assert topomesh.nb_regions(0,pid) == 0
topomesh.remove_wisp(0,pid)
return True
except AssertionError:
logging.info("".join([" " for l in range(loglevel)])+"Impossible to remove vertex : wrong face definition")
return False
[docs]def topomesh_remove_interface_edge(topomesh,eid, verbose=False, debug=False, loglevel=0):
edge_fids = list(topomesh.regions(1,eid))
fid_to_keep = np.min(edge_fids)
edge_fid_cells = [list(topomesh.regions(2,fid)) for fid in edge_fids]
edge_cells = np.unique(edge_fid_cells)
assert np.all([set(list(cids)) == set(list(edge_cells)) for cids in edge_fid_cells])
fids_to_remove = np.sort(edge_fids)[1:]
topomesh.unlink(2,fid_to_keep,eid)
for fid in fids_to_remove:
face_eids = list(topomesh.borders(2,fid))
for eid_to_link in face_eids:
topomesh.unlink(2,fid,eid_to_link)
if eid_to_link != eid and not eid_to_link in topomesh.borders(2,fid_to_keep):
topomesh.link(2,fid_to_keep,eid_to_link)
topomesh.remove_wisp(2,fid)
edge_pids = list(topomesh.borders(1,eid))
for pid in edge_pids:
topomesh.unlink(1,eid,pid)
if topomesh.nb_regions(0,pid) == 0:
topomesh.remove_wisp(0,pid)
topomesh.remove_wisp(1,eid)
return True
[docs]def topomesh_remove_boundary_vertex(topomesh, pid, verbose=False, debug=False, loglevel=0):
try:
vertex_eids = list(topomesh.regions(0,pid))
eid_to_keep = np.min(vertex_eids)
vertex_eid_faces = [list(topomesh.regions(1,eid)) for eid in vertex_eids]
vertex_faces = np.unique(vertex_eid_faces)
assert len(vertex_eids) == 2
assert np.all([set(list(fids)) == set(list(vertex_faces)) for fids in vertex_eid_faces])
assert len(np.unique([list(topomesh.borders(1,eid)) for eid in vertex_eids])) == 3
eids_to_remove = np.sort(vertex_eids)[1:]
topomesh.unlink(1,eid_to_keep,pid)
for eid in eids_to_remove:
edge_fids = list(topomesh.regions(1,eid))
for fid in edge_fids:
topomesh.unlink(2,fid,eid)
for pid_to_link in topomesh.borders(1,eid):
topomesh.unlink(1,eid,pid_to_link)
if pid_to_link != pid:
topomesh.link(1,eid_to_keep,pid_to_link)
topomesh.remove_wisp(1,eid)
topomesh.remove_wisp(0,pid)
# if np.array([list(topomesh.borders(1,eid)) for eid in topomesh.wisps(1)]).ndim != 2:
# logging.info("".join([" " for l in range(loglevel)])+eid_to_keep, eids_to_remove
# logging.info("".join([" " for l in range(loglevel)])+np.array(list(topomesh.wisps(1)))[np.array(map(len,[list(topomesh.borders(1,eid)) for eid in topomesh.wisps(1)]))!=2]
# raw_input()
return True
except AssertionError:
logging.info("".join([" " for l in range(loglevel)])+"Impossible to remove vertex : wrong edge definition")
return False
