# -*- coding: utf-8 -*-
# -*- python -*-
#
# PropertyTopomesh
#
# Copyright 2014-2016 INRIA - CIRAD - INRA
#
# File author(s): Guillaume Cerutti <guillaume.cerutti@inria.fr>,
# Hadrien Oliveri <hadrien.oliveri@inria.fr>
#
# File contributor(s): Guillaume Cerutti <guillaume.cerutti@inria.fr>,
# Hadrien Oliveri <hadrien.oliveri@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 I/O
=========================
Implementing functions to read and write PropertyTopomesh objects in several standard
mesh file formats, including PLY, OBJ or MSH.
.. autosummary::
:toctree: generated/
save_ply_property_topomesh -- write an ascii PLY file contaning a PropertyTopomesh
read_ply_property_topomesh -- reconstruct a PropertyTopomesh object from a PLY file
"""
from copy import copy
from time import time
import os
import re
import sys
import pickle
import numpy as np
from scipy import ndimage as nd
from cellcomplex.property_topomesh import PropertyTopomesh
from cellcomplex.property_topomesh.utils.matching_tools import kd_tree_match
from cellcomplex.property_topomesh.analysis import compute_topomesh_property
from cellcomplex.utils import array_dict
from sys import path
[docs]def save_property_topomesh(topomesh, path, cells_to_save=None, properties_to_save=dict([(0,['barycenter']),(1,[]),(2,[]),(3,[])]),**kwargs):
if cells_to_save is None:
cells_to_save = np.array(list(topomesh.wisps(3)))
triangles_to_save = np.array(np.unique(np.concatenate([np.array(list(topomesh.borders(3,c))) for c in cells_to_save])),int)
edges_to_save = np.array(np.unique(np.concatenate([np.array(list(topomesh.borders(2,t))) for t in triangles_to_save])),int)
vertices_to_save = np.array(np.unique(np.concatenate([np.array(list(topomesh.borders(1,e))) for e in edges_to_save])),int)
original_pids = kwargs.get('original_pids',False)
original_eids = kwargs.get('original_eids',False)
original_fids = kwargs.get('original_fids',False)
original_cids = kwargs.get('original_cids',True)
sub_topomesh = PropertyTopomesh(3)
vertices_to_pids = {}
for v in vertices_to_save:
if original_pids:
pid = sub_topomesh.add_wisp(0,v)
else:
pid = sub_topomesh.add_wisp(0)
vertices_to_pids[v] = pid
edges_to_eids = {}
for e in edges_to_save:
if original_eids:
eid = sub_topomesh.add_wisp(1,e)
else:
eid = sub_topomesh.add_wisp(1)
edges_to_eids[e] = eid
for v in topomesh.borders(1,e):
sub_topomesh.link(1,eid,vertices_to_pids[v])
triangles_to_fids = {}
for t in triangles_to_save:
if original_fids:
fid = sub_topomesh.add_wisp(2,t)
else:
fid = sub_topomesh.add_wisp(2)
triangles_to_fids[t] = fid
for e in topomesh.borders(2,t):
sub_topomesh.link(2,fid,edges_to_eids[e])
cells_to_cids = {}
for c in cells_to_save:
if original_cids:
cid = sub_topomesh.add_wisp(3,c)
else:
cid = sub_topomesh.add_wisp(3,c)
cells_to_cids[c] = cid
for t in topomesh.borders(3,c):
sub_topomesh.link(3,cid,triangles_to_fids[t])
wisps_to_save = {}
wisps_to_save[0] = vertices_to_save
wisps_to_save[1] = edges_to_save
wisps_to_save[2] = triangles_to_save
wisps_to_save[3] = cells_to_save
wisps_to_wids = {}
wisps_to_wids[0] = vertices_to_pids
wisps_to_wids[1] = edges_to_eids
wisps_to_wids[2] = triangles_to_fids
wisps_to_wids[3] = cells_to_cids
if not 0 in properties_to_save.keys():
properties_to_save[0] = []
if not 'barycenter' in properties_to_save[0]:
properties_to_save[0].append('barycenter')
for degree in properties_to_save.keys():
for property_name in properties_to_save[degree]:
print("Property ",property_name,'(',degree,')')
if topomesh.has_wisp_property(property_name,degree=degree,is_computed=True):
wids_to_save = array_dict(wisps_to_wids[degree]).values(wisps_to_save[degree])
sub_topomesh.update_wisp_property(property_name,degree,array_dict(topomesh.wisp_property(property_name,degree).values(wisps_to_save[degree]),wids_to_save))
pickle.dump(sub_topomesh,open(path,"wb"))
[docs]def save_tissue_property_topomesh(topomesh,tissue_filename='vplants.zip'):
from openalea.celltissue import Tissue, TissueDB, topen, Config, ConfigItem, ConfigFormat
tissue = Tissue()
ptyp = tissue.add_type("point")
etyp = tissue.add_type("edge")
ftyp = tissue.add_type("face")
ctyp = tissue.add_type("cell")
mesh_id = tissue.add_relation("mesh",(ptyp,etyp,ftyp,ctyp))
mesh = tissue.relation(mesh_id)
cells = {}
for c in topomesh.wisps(3):
if len(list(topomesh.borders(3,c)))>0:
cid = tissue.add_element(ctyp)
mesh.add_wisp(3,cid)
cells[c] = cid
faces = {}
for t in topomesh.wisps(2):
fid = tissue.add_element(ftyp)
mesh.add_wisp(2,fid)
faces[t] = fid
edges = {}
for e in topomesh.wisps(1):
eid = tissue.add_element(etyp)
mesh.add_wisp(1,eid)
edges[e] = eid
points = {}
for v in topomesh.wisps(0):
pid = tissue.add_element(ptyp)
mesh.add_wisp(0,pid)
points[v] = pid
for e in topomesh.wisps(1):
for v in topomesh.borders(1,e):
mesh.link(1,edges[e],points[v])
for t in topomesh.wisps(2):
for e in topomesh.borders(2,t):
mesh.link(2,faces[t],edges[e])
for c in topomesh.wisps(3):
if len(list(topomesh.borders(3,c)))>0:
for t in topomesh.borders(3,c):
mesh.link(3,cells[c],faces[t])
positions = {}
for v in topomesh.wisps(0):
positions[points[v]] = tuple(topomesh.wisp_property('barycenter',0)[v])
cfg = ConfigFormat("config")
cfg.add_section("topology")
cfg.add_item(ConfigItem("cell",ctyp) )
cfg.add_item(ConfigItem("face",ftyp) )
cfg.add_item(ConfigItem("edge",etyp) )
cfg.add_item(ConfigItem("point",ptyp) )
cfg.add_item(ConfigItem("mesh_id",mesh_id))
tissue_db = TissueDB()
tissue_db.set_tissue(tissue)
tissue_db.set_property("position",positions)
tissue_db.set_description("position","position of points")
tissue_db.set_config("config",cfg.config())
tissue_db.write(tissue_filename)
[docs]def multidim_data_ply_string(data):
"""Convert a multidimensional data into a PLY-formatted string
Parameters
----------
data : :class:`numpy.ndarray`
Returns
-------
"""
data_string = ""
data = np.array([l for l in data])
for dim in range(data.ndim):
data_string += " "+str(data.shape[dim])
# print data," (",len(data),")"
# data_string = " "+str(len(data))
for d in data.ravel():
data_string += " "+str(d)
# data_string += (" "+str(d)) if data.ndim==1 else multidim_data_ply_string(d)
return data_string
[docs]def save_ply_cellcomplex_topomesh(topomesh,ply_filename,properties_to_save=dict([(0,[]),(1,['length']),(2,['area','epidermis']),(3,[])]),color_faces=False,colormap=None,position_name='barycenter',oriented=True):
"""Implementing the PLY standard defined at Saisnbury Computational Workshop 2015.
Parameters
----------
topomesh
ply_filename
properties_to_save
color_faces
colormap
position_name
oriented
Returns
-------
"""
start_time =time()
print("--> Saving .ply")
ply_file = open(ply_filename,'w+')
if oriented:
compute_topomesh_property(topomesh,'oriented_borders',2)
compute_topomesh_property(topomesh,'oriented_vertices',2)
compute_topomesh_property(topomesh,'oriented_borders',3)
property_types = {}
property_types['bool'] = "int"
property_types['int'] = "int"
property_types['int32'] = "int"
property_types['int64'] = "int"
property_types['float'] = "float"
property_types['float32'] = "float"
property_types['float64'] = "float"
property_types['object'] = "list"
ply_file.write("ply\n")
ply_file.write("format ascii 1.0\n")
ply_file.write("element vertex "+str(topomesh.nb_wisps(0))+"\n")
ply_file.write("property float x\n")
ply_file.write("property float y\n")
ply_file.write("property float z\n")
for property_name in properties_to_save[0]:
if property_name != position_name and topomesh.has_wisp_property(property_name,0,is_computed=True):
property_type = property_types[str(topomesh.wisp_property(property_name,0).values().dtype)]
ply_file.write("property "+property_type+" "+property_name+"\n")
ply_file.write("element face "+str(topomesh.nb_wisps(2))+"\n")
ply_file.write("property list int int vertex_index\n")
if color_faces:
ply_file.write("property uchar red\n")
ply_file.write("property uchar green\n")
ply_file.write("property uchar blue\n")
triangle_data = array_dict(np.array([(list(topomesh.regions(2,t))[0])%256 for t in topomesh.wisps(2)]),list(topomesh.wisps(2)))
for property_name in properties_to_save[2]:
if topomesh.has_wisp_property(property_name,2,is_computed=True):
property_type = property_types[str(topomesh.wisp_property(property_name,2).values().dtype)]
ply_file.write("property "+property_type+" "+property_name+"\n")
ply_file.write("element edge "+str(topomesh.nb_wisps(1))+"\n")
ply_file.write("property int source\n")
ply_file.write("property int target\n")
for property_name in properties_to_save[1]:
if topomesh.has_wisp_property(property_name,1,is_computed=True):
property_type = property_types[str(topomesh.wisp_property(property_name,1).values().dtype)]
ply_file.write("property "+property_type+" "+property_name+"\n")
ply_file.write("element volume "+str(topomesh.nb_wisps(3))+"\n")
ply_file.write("property list int int face_index\n")
ply_file.write("property int label\n")
for property_name in properties_to_save[3]:
if topomesh.has_wisp_property(property_name,3,is_computed=True):
property_type = property_types[str(topomesh.wisp_property(property_name,3).values().dtype)]
ply_file.write("property "+property_type+" "+property_name+"\n")
# ply_file.write("element global 1\n")
# ply_file.write("property list char type\n")
ply_file.write("end_header\n")
vertex_index = {}
for v,pid in enumerate(topomesh.wisps(0)):
ply_file.write(str(topomesh.wisp_property('barycenter',0)[pid][0])+" ")
ply_file.write(str(topomesh.wisp_property('barycenter',0)[pid][1])+" ")
ply_file.write(str(topomesh.wisp_property('barycenter',0)[pid][2])+" ")
for property_name in properties_to_save[0]:
if property_name != position_name and topomesh.has_wisp_property(property_name,0,is_computed=True):
property_type = property_types[str(topomesh.wisp_property(property_name,0).values().dtype)]
if property_type == 'int':
ply_file.write(str(int(topomesh.wisp_property(property_name,0)[pid]))+" ")
else:
ply_file.write(str(topomesh.wisp_property(property_name,0)[pid])+" ")
ply_file.write("\n")
vertex_index[pid] = v
face_index = {}
for t,fid in enumerate(topomesh.wisps(2)):
if oriented:
oriented_face_pids = topomesh.wisp_property('oriented_vertices',2)[fid]
else:
face_pids = np.array(list(topomesh.borders(2,fid,2)))
face_edges = np.array([list(topomesh.borders(1,eid)) for eid in topomesh.borders(2,fid)])
# print fid," : ",face_pids, face_edges
oriented_face_pids = [face_pids[0]]
while len(oriented_face_pids) < len(face_pids):
current_pid = oriented_face_pids[-1]
pid_edges = face_edges[np.where(face_edges==current_pid)[0]]
candidate_pids = set(list(np.unique(pid_edges))).difference({current_pid})
if len(oriented_face_pids)>1:
candidate_pids = candidate_pids.difference({oriented_face_pids[-2]})
if len(list(candidate_pids))==0:
candidate_pids = set(list(face_pids)).difference(list(oriented_face_pids))
oriented_face_pids += [list(candidate_pids)[0]]
# print fid," : ",oriented_face_pids
ply_file.write(str(len(list(topomesh.borders(2,fid,2))))+" ")
for pid in oriented_face_pids:
ply_file.write(str(vertex_index[pid])+" ")
if color_faces:
if colormap is None:
ply_file.write(str(triangle_data[fid])+" ")
ply_file.write(str(triangle_data[fid])+" ")
ply_file.write(str(triangle_data[fid])+" ")
else:
color = colormap._color_points.values()[int(triangle_data[fid])]
ply_file.write(str(int(255*color[0]))+" ")
ply_file.write(str(int(255*color[1]))+" ")
ply_file.write(str(int(255*color[2]))+" ")
for property_name in properties_to_save[2]:
if topomesh.has_wisp_property(property_name,2,is_computed=True):
property_type = property_types[str(topomesh.wisp_property(property_name,2).values().dtype)]
if property_type == 'int':
ply_file.write(str(int(topomesh.wisp_property(property_name,2)[fid]))+" ")
else:
ply_file.write(str(topomesh.wisp_property(property_name,2)[fid])+" ")
ply_file.write("\n")
face_index[fid] = t
edge_index = {}
for e,eid in enumerate(topomesh.wisps(1)):
ply_file.write(str(vertex_index[list(topomesh.borders(1,eid))[0]])+" ")
ply_file.write(str(vertex_index[list(topomesh.borders(1,eid))[1]])+" ")
for property_name in properties_to_save[1]:
if topomesh.has_wisp_property(property_name,1,is_computed=True):
ply_file.write(str(topomesh.wisp_property(property_name,1)[eid])+" ")
ply_file.write("\n")
edge_index[eid] = e
for c, cid in enumerate(topomesh.wisps(3)):
if oriented:
oriented_cell_fids = topomesh.wisp_property('oriented_borders',3)[cid][0]
oriented_cell_fid_orientations = topomesh.wisp_property('oriented_borders',3)[cid][1]
ply_file.write(str(len(oriented_cell_fids))+" ")
for fid,orientation in zip(oriented_cell_fids,oriented_cell_fid_orientations):
ply_file.write(str(orientation*(face_index[fid]+1))+" ")
else:
for fid in topomesh.borders(3,cid):
ply_file.write(str(face_index[fid]+1)+" ")
ply_file.write(str(cid)+" ")
for property_name in properties_to_save[3]:
if topomesh.has_wisp_property(property_name,3,is_computed=True):
ply_file.write(str(topomesh.wisp_property(property_name,3)[cid])+" ")
ply_file.write("\n")
# ply_file.write("12 cell_complex\n")
ply_file.flush()
ply_file.close()
end_time = time()
print("<-- Saving .ply [",end_time-start_time,"s]")
[docs]def save_ply_property_topomesh(topomesh,ply_filename,properties_to_save=dict([(0,[]),(1,['length']),(2,['area','epidermis']),(3,[])]), position_name='barycenter', verbose=True):
"""Write a PropertyTopomesh as a standard PLY file.
Save a PropertyTopomesh object into a ASCII-formatted PLY file with a
subset of element properties written into the file. The written file is
readable by other mesh processing tools, and contains all the necessary
information to reconstruct the PropertyTopomesh including all topological
relationships and element properties.
Parameters
----------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
The PropertyTopomesh instance to save
ply_filename : str
Path to the PLY file to write to
properties_to_save : dict
Dictionary containing ,for each degree, the list of properties to write in the file
position_name : str
Property name used to define the XYZ position of vertices (default : `barycenter`)
verbose : bool
Whether to print information while writing the file
"""
if verbose:
start_time =time()
print("--> Saving .ply")
ply_file = open(ply_filename,'w+')
def get_property_type(dtype):
property_types = {}
property_types['bool'] = "int"
property_types['int'] = "int"
property_types['int16'] = "int"
property_types['int32'] = "int"
property_types['int64'] = "int"
property_types['uint8'] = "int"
property_types['uint16'] = "int"
property_types['float'] = "float"
property_types['float32'] = "float"
property_types['float64'] = "float"
property_types['object'] = "list"
if dtype in property_types:
return property_types[dtype]
elif re.match("[<|][SU][0-9]*", dtype):
return "str"
else:
raise KeyError("Data type not understood :"+str(dtype))
ply_file.write("ply\n")
ply_file.write("format ascii 1.0\n")
ply_file.write("element vertex "+str(topomesh.nb_wisps(0))+"\n")
ply_file.write("property float x\n")
ply_file.write("property float y\n")
ply_file.write("property float z\n")
for property_name in properties_to_save[0]:
if property_name != position_name and topomesh.has_wisp_property(property_name,0,is_computed=True):
property_data = topomesh.wisp_property(property_name,0).values()
if property_data.ndim == 1:
property_type = get_property_type(str(property_data.dtype))
if property_type == "list":
property_type = "list int "+get_property_type(str(np.array(property_data[0]).dtype))
elif property_data.ndim == 2:
property_type = "list int "+get_property_type(str(property_data.dtype))
if get_property_type(str(np.array(property_data).dtype)) == "list":
property_type = "tensor int int "+get_property_type(str(np.array(property_data[0][0]).dtype))
else:
property_type = "tensor "
for i in range(property_data.ndim - 1):
property_type += "int "
property_type += get_property_type(str(property_data.dtype))
ply_file.write("property "+property_type+" "+property_name+"\n")
ply_file.write("element face "+str(topomesh.nb_wisps(2))+"\n")
ply_file.write("property list int int vertex_index\n")
# if color_faces:
# ply_file.write("property uchar red\n")
# ply_file.write("property uchar green\n")
# ply_file.write("property uchar blue\n")
# triangle_data = array_dict(np.array([list(topomesh.regions(2,t))[0]%256 for t in topomesh.wisps(2)]),list(topomesh.wisps(2)))
for property_name in properties_to_save[2]:
if topomesh.has_wisp_property(property_name,2,is_computed=True):
property_data = topomesh.wisp_property(property_name,2).values()
if property_data.ndim == 1:
property_type = get_property_type(str(property_data.dtype))
if property_type == "list":
property_type = "list int "+get_property_type(str(np.array(property_data[0]).dtype))
elif property_data.ndim == 2:
property_type = "list int "+get_property_type(str(property_data.dtype))
if get_property_type(str(np.array(property_data).dtype)) == "list":
property_type = "tensor int int "+get_property_type(str(np.array(property_data[0][0]).dtype))
else:
property_type = "tensor "
for i in range(property_data.ndim - 1):
property_type += "int "
property_type += get_property_type(str(property_data.dtype))
ply_file.write("property "+property_type+" "+property_name+"\n")
ply_file.write("element edge "+str(topomesh.nb_wisps(1))+"\n")
ply_file.write("property int source\n")
ply_file.write("property int target\n")
ply_file.write("property list int int face_index\n")
for property_name in properties_to_save[1]:
if topomesh.has_wisp_property(property_name,1,is_computed=True):
property_data = topomesh.wisp_property(property_name,1).values()
if property_data.ndim == 1:
property_type = get_property_type(str(property_data.dtype))
if property_type == "list":
property_type = "list int "+get_property_type(str(np.array(property_data[0]).dtype))
elif property_data.ndim == 2:
property_type = "list int "+get_property_type(str(property_data.dtype))
if get_property_type(str(np.array(property_data).dtype)) == "list":
property_type = "tensor int int "+get_property_type(str(np.array(property_data[0][0]).dtype))
else:
property_type = "tensor "
for i in range(property_data.ndim - 1):
property_type += "int "
property_type += get_property_type(str(property_data.dtype))
ply_file.write("property "+property_type+" "+property_name+"\n")
ply_file.write("element volume "+str(topomesh.nb_wisps(3))+"\n")
ply_file.write("property list int int face_index\n")
ply_file.write("property int label\n")
for property_name in properties_to_save[3]:
if topomesh.has_wisp_property(property_name,3,is_computed=True):
property_data = topomesh.wisp_property(property_name,3).values()
if property_data.ndim == 1:
property_type = get_property_type(str(property_data.dtype))
if property_type == "list":
property_type = "list int "+get_property_type(str(np.array(property_data[0]).dtype))
elif property_data.ndim == 2:
property_type = "list int "+get_property_type(str(property_data.dtype))
if get_property_type(str(np.array(property_data).dtype)) == "list":
property_type = "tensor int int "+get_property_type(str(np.array(property_data[0][0]).dtype))
else:
property_type = "tensor "
for i in range(property_data.ndim - 1):
property_type += "int "
property_type += get_property_type(str(property_data.dtype))
ply_file.write("property "+property_type+" "+property_name+"\n")
ply_file.write("end_header\n")
vertex_index = {}
for v,pid in enumerate(topomesh.wisps(0)):
ply_file.write(str(topomesh.wisp_property(position_name,0)[pid][0])+" ")
ply_file.write(str(topomesh.wisp_property(position_name,0)[pid][1])+" ")
ply_file.write(str(topomesh.wisp_property(position_name,0)[pid][2]))
for property_name in properties_to_save[0]:
if property_name != position_name and topomesh.has_wisp_property(property_name,0,is_computed=True):
data = np.array(topomesh.wisp_property(property_name,0)[pid])
if data.ndim == 0:
ply_file.write(" "+str(data))
else:
ply_file.write(multidim_data_ply_string(data))
# for property_name in properties_to_save[0] :
# if property_name != position_name and topomesh.has_wisp_property(property_name,0,is_computed=True):
# property_type = property_types[str(topomesh.wisp_property(property_name,0).values().dtype)]
# if topomesh.wisp_property(property_name,0).values().ndim == 1:
# if property_type == 'int':
# ply_file.write(str(int(topomesh.wisp_property(property_name,0)[pid]))+" ")
# else:
# ply_file.write(str(topomesh.wisp_property(property_name,0)[pid])+" ")
# elif topomesh.wisp_property(property_name,0).values().ndim == 2:
# ply_file.write(str(int(len(topomesh.wisp_property(property_name,0)[fid])))+" ")
# for p in topomesh.wisp_property(property_name,0)[fid]:
# if property_type == 'int':
# ply_file.write(str(int(p))+" ")
# else:
# ply_file.write(str(p)+" ")
ply_file.write("\n")
vertex_index[pid] = v
face_index = {}
for t,fid in enumerate(topomesh.wisps(2)):
ply_file.write(str(len(list(topomesh.borders(2,fid,2))))+" ")
if topomesh.has_wisp_property('oriented_vertices',2,is_computed=True):
for i_pid, pid in enumerate(topomesh.wisp_property('oriented_vertices',2)[fid]):
ply_file.write((" " if i_pid>0 else "")+str(vertex_index[pid]))
else:
for i_pid, pid in enumerate(topomesh.borders(2,fid,2)):
ply_file.write((" " if i_pid>0 else "")+str(vertex_index[pid]))
# ply_file.write(str(len(list(topomesh.borders(2,fid))))+" ")
# for eid in topomesh.borders(2,fid):
# ply_file.write(str(edge_index[eid])+" ")
# if color_faces:
# ply_file.write(str(triangle_data[fid])+" ")
# ply_file.write(str(triangle_data[fid])+" ")
# ply_file.write(str(triangle_data[fid])+" ")
# if colormap is None:
# ply_file.write(str(triangle_data[fid])+" ")
# ply_file.write(str(triangle_data[fid])+" ")
# ply_file.write(str(triangle_data[fid])+" ")
# else:
# color = colormap._color_points.values()[int(triangle_data[fid])]
# ply_file.write(str(int(255*color[0]))+" ")
# ply_file.write(str(int(255*color[1]))+" ")
# ply_file.write(str(int(255*color[2]))+" ")
for property_name in properties_to_save[2]:
if topomesh.has_wisp_property(property_name, 2, is_computed=True):
data = np.array(topomesh.wisp_property(property_name, 2)[fid])
if data.ndim == 0:
ply_file.write(" " + str(data))
else:
ply_file.write(multidim_data_ply_string(data))
# for property_name in properties_to_save[2]:
# if topomesh.has_wisp_property(property_name,2,is_computed=True):
# property_type = property_types[str(topomesh.wisp_property(property_name,2).values().dtype)]
# if topomesh.wisp_property(property_name,2).values().ndim == 1:
# if property_type == 'int':
# ply_file.write(str(int(topomesh.wisp_property(property_name,2)[fid]))+" ")
# else:
# ply_file.write(str(topomesh.wisp_property(property_name,2)[fid])+" ")
# elif topomesh.wisp_property(property_name,2).values().ndim == 2:
# ply_file.write(str(int(len(topomesh.wisp_property(property_name,2)[fid])))+" ")
# for p in topomesh.wisp_property(property_name,2)[fid]:
# if property_type == 'int':
# ply_file.write(str(int(p))+" ")
# else:
# ply_file.write(str(p)+" ")
ply_file.write("\n")
face_index[fid] = t
edge_index = {}
for e,eid in enumerate(topomesh.wisps(1)):
ply_file.write(str(vertex_index[list(topomesh.borders(1,eid))[0]])+" ")
ply_file.write(str(vertex_index[list(topomesh.borders(1,eid))[1]])+" ")
ply_file.write(str(len(list(topomesh.regions(1,eid))))+" ")
for i_fid, fid in enumerate(topomesh.regions(1,eid)):
ply_file.write((" " if i_fid>0 else "")+str(face_index[fid]))
for property_name in properties_to_save[1]:
if topomesh.has_wisp_property(property_name, 1, is_computed=True):
data = np.array(topomesh.wisp_property(property_name, 1)[eid])
if data.ndim == 0:
ply_file.write(" " + str(data))
else:
ply_file.write(multidim_data_ply_string(data))
# for property_name in properties_to_save[1]:
# if topomesh.has_wisp_property(property_name,1,is_computed=True):
# property_type = property_types[str(topomesh.wisp_property(property_name,1).values().dtype)]
# if topomesh.wisp_property(property_name,1).values().ndim == 1:
# if property_type == 'int':
# ply_file.write(str(int(topomesh.wisp_property(property_name,1)[eid]))+" ")
# else:
# ply_file.write(str(topomesh.wisp_property(property_name,1)[eid])+" ")
# elif topomesh.wisp_property(property_name,1).values().ndim == 2:
# ply_file.write(str(int(len(topomesh.wisp_property(property_name,1)[eid])))+" ")
# for p in topomesh.wisp_property(property_name,1)[eid]:
# if property_type == 'int':
# ply_file.write(str(int(p))+" ")
# else:
# ply_file.write(str(p)+" ")
ply_file.write("\n")
edge_index[eid] = e
for c, cid in enumerate(topomesh.wisps(3)):
ply_file.write(str(len(list(topomesh.borders(3,cid))))+" ")
for i_fid, fid in enumerate(topomesh.borders(3,cid)):
ply_file.write((" " if i_fid>0 else "")+str(face_index[fid]))
ply_file.write(" "+str(cid))
for property_name in properties_to_save[3]:
if topomesh.has_wisp_property(property_name, 3, is_computed=True):
data = np.array(topomesh.wisp_property(property_name, 3)[cid])
if data.ndim == 0:
ply_file.write(" " + str(data))
else:
ply_file.write(multidim_data_ply_string(data))
# if topomesh.has_wisp_property(property_name,3,is_computed=True):
# property_type = property_types[str(topomesh.wisp_property(property_name,3).values().dtype)]
# if topomesh.wisp_property(property_name,3).values().ndim == 1:
# if property_type == 'int':
# ply_file.write(str(int(topomesh.wisp_property(property_name,3)[cid]))+" ")
# else:
# ply_file.write(str(topomesh.wisp_property(property_name,3)[cid])+" ")
# elif topomesh.wisp_property(property_name,3).values().ndim == 2:
# ply_file.write(str(int(len(topomesh.wisp_property(property_name,3)[cid])))+" ")
# for p in topomesh.wisp_property(property_name,3)[cid]:
# if property_type == 'int':
# ply_file.write(str(int(p))+" ")
# else:
# ply_file.write(str(p)+" ")
ply_file.write("\n")
ply_file.flush()
ply_file.close()
if verbose:
end_time = time()
print("<-- Saving .ply [",end_time-start_time,"s]")
[docs]def read_ply_property_topomesh(ply_filename, verbose=False, fuse_elements=True, fuse_vertices=True):
"""Read a PropertyTopomesh stored as a PLY file.
Read an ASCII-encoded PLY file containing a cell complex. The file can be a
regular mesh file with only `vertex` and `face` elements, but it can also
contain `edge` and `volume` elements that define topological elements of
dimensions 1 and 3 respectively. Associated properties can be of int, float,
string, vector or tensor types.
Parameters
----------
ply_filename : str
Path to the PLY file to read from
verbose : bool
Whether to print all information read from the file
fuse_elements : bool
Whether identical topological elements should be fused in a single wisp
fuse_vertices : bool
Whether vertices at the same position should be fused in a single wisp
Returns
-------
topomesh : :class:`cellcomplex.property_topomesh.PropertyTopomesh`
Cell complex with all the element properties read from the file
"""
if verbose:
start_time =time()
print("--> Reading .ply")
property_types = {}
property_types['int'] = 'int'
property_types['int32'] = 'int'
property_types['int16'] = 'int'
property_types['uint'] = 'int'
property_types['uint8'] = 'int'
property_types['uint16'] = 'int'
property_types['uchar'] = 'int'
property_types['float'] = 'float'
property_types['float32'] = 'float'
property_types['list'] = 'list'
property_types['tensor'] = 'tensor'
property_types['string'] = 'str'
property_types['|S7'] = 'str'
property_types['str'] = 'str'
ply_file = open(ply_filename,'rU')
ply_stream = enumerate(ply_file,1)
assert "ply" in next(ply_stream)[1]
assert "ascii" in next(ply_stream)[1]
n_wisps = {}
properties = {}
properties_types = {}
properties_list_types = {}
properties_tensor_dims = {}
element_name = ""
property_name = ""
elements = []
lineno, line = next(ply_stream)
while not 'end_header' in line:
try:
if re.split(' ',line)[0] == 'element':
element_name = re.split(' ',line)[1]
elements.append(element_name)
n_wisps[element_name] = int(re.split(' ',line)[2])
properties[element_name] = []
properties_types[element_name] = {}
properties_list_types[element_name] = {}
properties_tensor_dims[element_name] = {}
if re.split(' ',line)[0] == 'property':
# print line
property_name = re.split(' ',line)[-1][:-1]
properties[element_name].append(property_name)
properties_types[element_name][property_name] = re.split(' ',line)[1]
if properties_types[element_name][property_name] == 'list':
list_type = re.split(' ',line)[-2]
properties_list_types[element_name][property_name] = list_type
elif properties_types[element_name][property_name] == 'tensor':
properties_tensor_dims[element_name][property_name] = (np.array(re.split(' ',line))[:-2]=='int').sum()
list_type = re.split(' ',line)[-2]
properties_list_types[element_name][property_name] = list_type
except Exception as e:
raise ValueError(ply_filename, lineno, line, e)
lineno, line = next(ply_stream)
if verbose: print(n_wisps)
if verbose: print(properties)
if verbose: print(properties_list_types)
element_properties = {}
for element_name in elements:
element_properties[element_name] = {}
for wid in range(n_wisps[element_name]):
lineno, line = next(ply_stream)
line = line.strip()
line_props = {}
prop_index = 0
try:
for prop in properties[element_name]:
prop_type = properties_types[element_name][prop]
if verbose:
print(element_name," ",wid," - ",prop," (",property_types[prop_type],")")
if property_types[prop_type] == 'float':
line_props[prop] = float(re.split(' ',line)[prop_index])
prop_index += 1
elif property_types[prop_type] == 'int':
line_props[prop] = int(re.split(' ',line)[prop_index])
prop_index += 1
elif property_types[prop_type] == 'str':
line_props[prop] = str(re.split(' ',line)[prop_index])
prop_index += 1
elif property_types[prop_type] == 'list':
list_length = int(re.split(' ',line)[prop_index])
prop_index += 1
list_type = properties_list_types[element_name][prop]
if property_types[list_type] == 'float':
line_props[prop] = [float(p) for p in re.split(' ',line)[prop_index:prop_index+list_length]]
elif property_types[list_type] == 'int':
line_props[prop] = [int(p) for p in re.split(' ',line)[prop_index:prop_index+list_length]]
prop_index += list_length
elif property_types[prop_type] == 'tensor':
n_dims = properties_tensor_dims[element_name][prop]
tensor_dims = tuple(np.array(re.split(' ',line)[prop_index:prop_index+n_dims]).astype(int))
prop_index += n_dims
list_type = properties_list_types[element_name][prop]
line_props[prop] = np.array(re.split(' ',line)[prop_index:prop_index+np.prod(tensor_dims)]).astype(property_types[list_type]).reshape(tensor_dims)
prop_index += np.prod(tensor_dims)
if verbose:
print(element_name," ",wid," - ",prop," : ",line_props[prop])
except Exception as e:
raise ValueError(ply_filename, lineno, line, prop, property_types[prop_type], e)
element_properties[element_name][wid] = line_props
ply_file.close()
if verbose:
print("<-- Parsing .ply [",time()-start_time,"s]")
element_matching = {}
point_positions = {}
for pid in range(n_wisps['vertex']):
point_positions[pid] = np.array([element_properties['vertex'][pid][dim] for dim in ['x','y','z']])
face_vertices = {}
for fid in range(n_wisps['face']):
if 'vertex_index' in element_properties['face'][fid]:
face_vertices[fid] = element_properties['face'][fid]['vertex_index']
elif 'vertex_indices' in element_properties['face'][fid]:
face_vertices[fid] = element_properties['face'][fid]['vertex_indices']
timecheck = verbose
if fuse_elements:
forcetest = False
if timecheck: check_time =time()
if verbose: print(len(point_positions)," Points, ", len(face_vertices), " Faces")
if fuse_vertices:
unique_points = np.unique(np.array(list(point_positions.values())),axis=0)
if not forcetest and len(unique_points) == len(point_positions):
unique_points = np.array(list(point_positions.values()))
point_matching = array_dict(point_positions.keys(),point_positions.keys())
else:
point_matching = array_dict(kd_tree_match(np.array(list(point_positions.values())),unique_points),list(point_positions.keys()))
else:
unique_points = np.array(list(point_positions.values()))
point_matching = array_dict(point_positions.keys(), point_positions.keys())
element_matching['vertex'] = point_matching
if verbose: print(len(unique_points)," Unique Points")
if timecheck: print('point unicity test:',time()-check_time)
if timecheck: check_time =time()
faces = np.array(list(face_vertices.values()))
if len(faces)==0:
if verbose: print("No faces !")
triangular = True
unique_triangles = []
triangle_matching = {}
elif faces.ndim == 2:
if verbose: print("Faces = Triangles !")
triangular = True
triangles = np.sort(point_matching.values(faces))
unique_triangles = np.unique(triangles,axis=0)
if not forcetest and len(unique_triangles) == len(triangles):
unique_triangles = triangles
triangle_matching = array_dict(face_vertices.keys(),face_vertices.keys())
else:
triangle_matching = array_dict(kd_tree_match(triangles,unique_triangles),face_vertices.keys())
else:
triangular = False
if len(faces)>0:
unique_triangles = point_matching.values(faces)
triangle_matching = array_dict(list(face_vertices.keys()),list(face_vertices.keys()))
else:
unique_triangles = np.array([])
triangle_matching = array_dict()
element_matching['face'] = triangle_matching
if verbose: print(len(unique_triangles)," Unique Faces")
if timecheck: print('face unicity test:',time()-check_time)
if timecheck: check_time =time()
if 'edge' in n_wisps:
edge_vertices = {}
edge_faces = {}
for eid in range(n_wisps['edge']):
edge_faces[eid] = []
edge_vertices[eid] = point_matching.values(np.array([element_properties['edge'][eid][dim] for dim in ['source','target']]))
if 'face_index' in element_properties['edge'][eid]:
edge_faces[eid] = element_properties['edge'][eid]['face_index']
#print element_properties['edge']
if not 'face_index' in properties['edge'] and len(faces)>0:
face_edge_vertices = np.sort(np.concatenate([np.transpose([v,list(v[1:])+[v[0]]]) for v in unique_triangles]))
face_edge_faces = np.concatenate([fid*np.ones_like(unique_triangles[fid]) for fid in range(len(unique_triangles))])
face_edge_matching = kd_tree_match(face_edge_vertices,np.sort(list(edge_vertices.values())))
for eid in range(n_wisps['edge']):
edge_faces[eid] = []
for eid, fid in zip(face_edge_matching, face_edge_faces):
if eid is not None:
edge_faces[eid] += [fid]
else:
if triangular:
if len(unique_triangles) == 0:
edge_vertices = {}
edge_faces = {}
else:
edge_vertices = dict(zip(list(range(3*len(unique_triangles))),np.sort(np.concatenate([np.transpose([v,list(v[1:])+[v[0]]]) for v in unique_triangles]))))
edge_faces = dict(zip(list(range(3*len(unique_triangles))),np.concatenate([fid*np.ones_like(unique_triangles[fid]) for fid in range(len(unique_triangles))])))
else:
edge_index = list(range(len(unique_triangles[0])))
for v in unique_triangles[1:]:
edge_index += [edge_index[-1]+1+r for r in range(len(v))]
edge_vertices = dict(list(zip(edge_index,np.sort(np.concatenate([np.transpose([v,list(v[1:])+[v[0]]]) for v in unique_triangles])))))
edge_faces = dict(list(zip(edge_index,np.concatenate([fid*np.ones_like(unique_triangles[fid]) for fid in range(len(unique_triangles))]))))
if verbose: print(len(edge_vertices)," Edges")
if len(edge_vertices)>0:
unique_edges = np.unique(np.sort(list(edge_vertices.values())),axis=0)
if not forcetest and len(unique_edges) == len(edge_vertices.values()):
unique_edges = list(edge_vertices.values())
edge_matching = array_dict(list(edge_vertices.keys()),list(edge_vertices.keys()))
else:
# edge_matching = array_dict(vq(np.sort(edge_vertices.values()),unique_edges)[0],edge_vertices.keys())
edge_matching = array_dict(kd_tree_match(np.sort(list(edge_vertices.values())),unique_edges),list(edge_vertices.keys()))
else:
unique_edges = np.array([])
edge_matching = array_dict()
element_matching['edge'] = edge_matching
if verbose: print(len(unique_edges)," Unique Edges")
if timecheck: print('edge unicity test:',time()-check_time)
if timecheck: check_time =time()
face_cells = {}
for fid in range(len(unique_triangles)):
face_cells[fid] = set()
cell_matching = {}
if 'volume' in n_wisps:
oriented_cells = False
if np.any([np.any([f<0 for f in element_properties['volume'][c]['face_index']]) for c in range(n_wisps['volume'])]):
oriented_cells = True
if 'label' in properties['volume']:
for cid in range(n_wisps['volume']):
cell_matching[cid] = element_properties['volume'][cid]['label']
else:
cell_matching = dict(list(zip(list(range(n_wisps['volume'])),list(range(n_wisps['volume'])))))
if 'face_index' in properties['volume']:
for c in range(n_wisps['volume']):
for f in element_properties['volume'][c]['face_index']:
if oriented_cells:
face_cells[triangle_matching[np.abs(f)-1]] = face_cells[triangle_matching[np.abs(f)-1]].union({cell_matching[c]})
else:
face_cells[triangle_matching[f]] = face_cells[triangle_matching[f]].union({cell_matching[c]})
else:
for f in range(len(unique_triangles)):
face_cells[triangle_matching[f]] = {0}
else:
cell_matching[0] = 0
for f in range(len(unique_triangles)):
face_cells[triangle_matching[f]] = {0}
element_matching['volume'] = cell_matching
if verbose: print(len(cell_matching)," Cells")
if timecheck: print('cell unicity test:',time()-check_time)
else:
unique_points = np.array(list(point_positions.values()))
point_matching = array_dict(list(point_positions.keys()),list(point_positions.keys()))
faces = np.array(face_vertices.values())
if len(faces) == 0:
if verbose: print("No faces !")
triangular = True
unique_triangles = []
triangle_matching = {}
elif faces.ndim == 2:
if verbose: print("Faces = Triangles !")
triangular = True
unique_triangles = np.sort(faces)
triangle_matching = array_dict(list(face_vertices.keys()),list(face_vertices.keys()))
else:
triangular = False
if len(faces)>0:
unique_triangles = point_matching.values(faces)
triangle_matching = array_dict(list(face_vertices.keys()),list(face_vertices.keys()))
else:
unique_triangles = np.array([])
triangle_matching = array_dict()
element_matching['face'] = triangle_matching
if triangular:
if len(unique_triangles) == 0:
edge_vertices = {}
edge_faces = {}
else:
edge_vertices = dict(zip(list(range(3*len(unique_triangles))),np.sort(np.concatenate([np.transpose([v,list(v[1:])+[v[0]]]) for v in unique_triangles]))))
edge_faces = dict(zip(list(range(3*len(unique_triangles))),np.concatenate([fid*np.ones_like(unique_triangles[fid]) for fid in range(len(unique_triangles))])))
else:
edge_index = list(range(len(unique_triangles[0])))
for v in unique_triangles[1:]:
edge_index += [edge_index[-1]+1+r for r in range(len(v))]
edge_vertices = dict(list(zip(edge_index,np.sort(np.concatenate([np.transpose([v,list(v[1:])+[v[0]]]) for v in unique_triangles])))))
edge_faces = dict(list(zip(edge_index,np.concatenate([fid*np.ones_like(unique_triangles[fid]) for fid in range(len(unique_triangles))]))))
unique_edges = edge_vertices.values()
if len(unique_edges) == 0:
edge_matching = array_dict()
else:
edge_matching = array_dict(edge_vertices.keys(),edge_vertices.keys())
element_matching['edge'] = edge_matching
face_cells = {}
for fid in range(len(unique_triangles)):
face_cells[fid] = set()
cell_matching = {}
if 'volume' in n_wisps:
if 'label' in properties['volume']:
for cid in range(n_wisps['volume']):
cell_matching[cid] = element_properties['volume'][cid]['label']
else:
cell_matching = dict(list(zip(list(range(n_wisps['volume'])),list(range(n_wisps['volume'])))))
if 'face_index' in properties['volume']:
for c in range(n_wisps['volume']):
for f in element_properties['volume'][c]['face_index']:
face_cells[triangle_matching[f]] = face_cells[triangle_matching[f]].union({cell_matching[c]})
else:
for f in range(len(unique_triangles)):
face_cells[triangle_matching[f]] = {0}
else:
cell_matching[0] = 0
for f in range(len(unique_triangles)):
face_cells[triangle_matching[f]] = {0}
element_matching['volume'] = cell_matching
if timecheck: check_time =time()
topomesh = PropertyTopomesh(3)
for pid in range(len(unique_points)):
topomesh.add_wisp(0,pid)
for fid in range(len(unique_triangles)):
topomesh.add_wisp(2,fid)
for cid in face_cells[fid]:
if not topomesh.has_wisp(3, cid):
topomesh.add_wisp(3,cid)
topomesh.link(3,cid,fid)
for eid,e in enumerate(unique_edges):
topomesh.add_wisp(1,eid)
for pid in e:
topomesh.link(1,eid,pid)
for e in edge_vertices.keys():
eid = edge_matching[e]
if isinstance(edge_faces[e],list):
for f in edge_faces[e]:
fid = triangle_matching[f]
topomesh.link(2,fid,eid)
else:
fid = edge_faces[e]
topomesh.link(2,fid,eid)
if timecheck: print('topomesh topo creation:',time()-check_time)
topomesh.update_wisp_property('barycenter',0,array_dict(unique_points,np.arange(len(unique_points))))
native_properties = {}
native_properties['vertex'] = ['x','y','z']
native_properties['edge'] = ['source','target','face_index']
native_properties['face'] = ['vertex_index','vertex_indices','red','green','blue']
native_properties['volume'] = ['face_index','label']
for property_degree, element_name in zip([0,1,2,3],['vertex','edge','face','volume']):
if element_name in n_wisps:
for property_name in properties[element_name]:
if property_name not in native_properties[element_name]:
if verbose: print("Adding "+property_name+" "+element_name+" property")
property_dict = {}
for w in element_properties[element_name].keys():
property_dict[element_matching[element_name][w]] = element_properties[element_name][w][property_name]
try:
topomesh.update_wisp_property(property_name,property_degree,array_dict(property_dict))
except ValueError:
property_dict = dict(zip(property_dict.keys(),[list(property_dict[k]) for k in property_dict.keys()]))
topomesh.update_wisp_property(property_name,property_degree,array_dict(property_dict))
if timecheck: print('topomesh creation:',time()-check_time)
if verbose:
end_time = time()
print("<-- Reading .ply [",end_time-start_time,"s]")
return topomesh
[docs]def read_obj_property_topomesh(obj_filename, verbose=True):
"""
"""
import re
# from cellcomplex.property_topomesh.utils.array_tools import array_unique
from cellcomplex.property_topomesh.creation import triangle_topomesh
if verbose:
start_time =time()
print("--> Reading .obj")
obj_file = open(obj_filename,'rU')
if not obj_file:
raise IOError("Unable to open "+obj_filename+"!")
obj_stream = enumerate(obj_file,1)
vertices = []
faces = []
vertex_normals = []
face_normals = []
for lineno, line in obj_stream:
if len(line)>0:
if re.split(' ',line)[0] == 'v':
vertices += [np.array(re.split(' ',line)[1:4]).astype(float)]
if re.split(' ',line)[0] == 'f':
if '//' in line:
faces += [np.array([re.split('//',v)[0] for v in re.split(' ',line)[1:]]).astype(int) - 1]
elif '/' in line:
faces += [np.array([re.split('/',v)[0] for v in re.split(' ',line)[1:]]).astype(int) - 1]
else:
faces += [np.array([v for v in re.split(' ',line)[1:4]]).astype(int) - 1]
vertex_positions = dict(list(zip(list(range(len(vertices))),vertices)))
topomesh = triangle_topomesh(faces,vertex_positions)
return topomesh
[docs]def save_obj_property_topomesh(topomesh, obj_filename, reorient_faces=False, verbose=False):
if reorient_faces:
compute_topomesh_property(topomesh,'normal',2,normal_method='orientation')
start_time =time()
print("--> Saving .obj + .mtl")
mtlfile = open(obj_filename[:-4]+".mtl",'w')
for c in topomesh.wisps(3):
mtlfile.write("newmtl mat"+str(c)+"\n")
col = [np.random.rand(),np.random.rand(),np.random.rand()]
mtlfile.write(" Ka "+str(col[0])+" "+str(col[1])+" "+str(col[2])+"\n")
mtlfile.write(" Kd "+str(col[0])+" "+str(col[1])+" "+str(col[2])+"\n")
mtlfile.write(" Ks 0.0 0.0 0.0\n")
mtlfile.write(" Ns 0.0\n")
mtlfile.write(" d 1.0\n")
mtlfile.write(" Tr 1.0\n")
objfile = open(obj_filename,'w')
objfile.write("#\n")
objfile.write("# Wavefront OBJ file\n")
objfile.write("# TextureMontage Atlas\n")
objfile.write("#\n")
objfile.write("mtllib "+obj_filename[:-4]+".mtl\n")
vertex_index = {}
for i,v in enumerate(list(topomesh.wisps(0))):
pos = topomesh.wisp_property('barycenter',0)[v]
objfile.write("v "+str(pos[0])+" "+str(pos[1])+" "+str(pos[2])+"\n")
vertex_index[v] = i
for c in topomesh.wisps(3):
objfile.write("g object"+str(c)+"\n")
objfile.write("usemtl mat"+str(c)+"\n")
for t in topomesh.borders(3,c):
if reorient_faces:
# tri = list(topomesh.borders(2,t,2))[::topomesh.wisp_property('orientation',2)[t]]
tri = list(topomesh.wisp_property('oriented_vertices',2)[t])[::topomesh.wisp_property('orientation',2)[t]]
else:
tri = list(topomesh.borders(2,t,2))
objfile.write("f "+str(vertex_index[tri[0]]+1)+"/"+str(c)+" "+str(vertex_index[tri[1]]+1)+"/"+str(c)+" "+str(vertex_index[tri[2]]+1)+"/"+str(c)+"\n")
objfile.close()
end_time = time()
print("<-- Saving .obj + .mtl [",end_time-start_time,"s]")
[docs]def read_msh_property_topomesh(msh_filename, preserve_cells=True, verbose=False):
"""
"""
import re
# from cellcomplex.property_topomesh.utils.array_tools import array_unique
from cellcomplex.property_topomesh.creation import triangle_topomesh
if verbose:
start_time =time()
print("--> Reading .msh")
msh_file = open(msh_filename,'rU')
if not msh_file:
raise IOError("Unable to open "+msh_filename+"!")
msh_stream = enumerate(msh_file,1)
vertex_positions = {}
faces = []
face_cells = []
vertex_normals = []
face_normals = []
lineno, line = next(msh_stream)
while not "$MeshFormat" in line:
lineno, line = next(msh_stream)
lineno, line = next(msh_stream)
mesh_format = line
while not "$Nodes" in line:
lineno, line = next(msh_stream)
lineno, line = next(msh_stream)
n_nodes = int(re.split(' ',line)[0])
print(n_nodes," Nodes")
for n in range(n_nodes):
lineno, line = next(msh_stream)
nid = int(re.split(' ',line)[0])
n_pos = np.array(re.split(' ',line)[1:4]).astype(float)
vertex_positions[nid] = n_pos
lineno, line = next(msh_stream)
assert "$EndNodes" in line
while not "$Elements" in line:
lineno, line = next(msh_stream)
lineno, line = next(msh_stream)
n_elements = int(re.split(' ',line)[0])
print(n_elements," Elements")
# 2 : triangular face
# 1 : edge
# 15 : point
element_node_number = dict([(2,3),(1,2),(15,1)])
#element_node_number = dict([(1,2),(15,1)])
for e in range(n_elements):
lineno, line = next(msh_stream)
eid = int(re.split(' ',line)[0])
e_type = int(re.split(' ',line)[1])
n_tags = int(re.split(' ',line)[2])
e_tags = np.array(re.split(' ',line)[3:3+n_tags]).astype(int)
cell = e_tags[-1]
try:
e_nodes = np.array([re.split(' ',line)[3+n_tags+i] for i in range(element_node_number[e_type])]).astype(int)
except KeyError:
raise KeyError("MSH Element type "+str(e_type)+" not handled yet!")
if e_type == 2:
faces += [list(e_nodes)]
face_cells += [cell]
lineno, line = next(msh_stream)
assert "$EndElements" in line
msh_file.close()
topomesh = triangle_topomesh(faces,vertex_positions)
if preserve_cells:
topomesh.remove_wisp(3,0)
for cid in np.sort(np.unique(face_cells)):
topomesh.add_wisp(3,cid)
for fid,cid in enumerate(face_cells):
topomesh.link(3,cid,fid)
return topomesh
[docs]def meshread(filename):
file_formats = ['.ply','.obj','.msh']
# extension = path(filename).ext
extension = os.path.splitext(filename)[1]
try:
assert extension in file_formats
except AssertionError:
raise IOError("The file format "+extension+" is not recognized, only "+str(file_formats)+" mesh files are supported")
else:
if extension == '.ply':
return read_ply_property_topomesh(filename)
elif extension == '.obj':
return read_obj_property_topomesh(filename)
elif extension == '.msh':
return read_msh_property_topomesh(filename)
