---
title: VISI601_CMI Maillages en Python et numpy
---
# Maillages en Python et numpy
> [name=Jacques-Olivier Lachaud][time=January 2022][color=#907bf7]
> (Les images peuvent être soumises à des droits d'auteur. Elles sont utilisées ici exclusivement dans un but pédagogique)
###### tags: `visi601`
Retour à [VISI601_CMI (Main) Algorithmique numérique](https://codimd.math.cnrs.fr/s/IWTaBkA9m)
---
[TOC]
---
Il y a beaucoup de façon de représenter des maillages (surfaces triangulées ou polygonales) en Python. Nous en proposons ici une qui a l'avantage d'être simple (en gros des tableaux `numpy`) et d'offrir une interface très pratique et évolutive [Polyscope - python](https://polyscope.run/py).
## Installation
On suppose que vous avez déjà installé [numpy](https://numpy.org/doc/stable/user/index.html). Pour installer Polyscope, rien de plus simple:
```python
python -m pip install polyscope
```
## Utilisation
Dans vos code importez `numpy`et `polyscope` avant d'utiliser les fonctions associées
```python
import polyscope as ps
import numpy as np
# Initialise polyscope
ps.init()
```
## Créer une surface avec des faces polygonales
Voici comment créer une pyramide avec une base carrée.
```python
# `verts` is a Nx3 numpy array of vertex positions
# `faces` is a Fx3 array of indices, or a nested list
verts=np.array([[1.,0.,0.],[0.,1.,0.],[-1.,0.,0.],[0.,-1.,0.],[0.,0.,1.]])
faces=[[0,1,2,3],[1,0,4],[2,1,4],[3,2,4],[0,3,4]]
ps.register_surface_mesh("my mesh", verts, faces )
```
## Donner la main à Polyscope pour visualiser la surface
Il suffit d'appeler la fonction `show()` de Polyscope et il affiche automatiquement toutes les surfaces qui ont été déclarées.
```python
ps.show()
```
## Charger une surface définie dans le format OBJ (Wavefront)
OBJ est un format de fichier relativement simple pour représenter des objets 3D. On va utiliser un code tout fait pour charger et/ou sauvegarder de telles données. Il est à 90% pris de http://jamesgregson.ca/loadsave-wavefront-obj-files-in-python.html, merci à lui. Copiez le code ci-dessous dans un fichier `wavefront.py`.
```python
# wavefront.py
import numpy as np
class WavefrontOBJ:
def __init__( self, default_mtl='default_mtl' ):
self.path = None # path of loaded object
self.mtllibs = [] # .mtl files references via mtllib
self.mtls = [ default_mtl ] # materials referenced
self.mtlid = [] # indices into self.mtls for each polygon
self.vertices = [] # vertices as an Nx3 or Nx6 array (per vtx colors)
self.normals = [] # normals
self.texcoords = [] # texture coordinates
self.polygons = [] # M*Nv*3 array, Nv=# of vertices, stored as\ vid,tid,nid (-1 for N/A)
def only_coordinates( self ):
V=np.zeros( ( len( self.vertices ), 3 ), np.float64 )
for i in range(len(self.vertices)):
V[i][0]= self.vertices[i][0]
V[i][1]= self.vertices[i][1]
V[i][2]= self.vertices[i][2]
return V
def only_faces( self ):
all_faces=[]
for f in self.polygons:
face=[]
for indices in f:
face.append( indices[0] )
all_faces.append( face )
return all_faces
def load_obj( filename: str, default_mtl='default_mtl', triangulate=False ) -> WavefrontOBJ:
"""Reads a .obj file from disk and returns a WavefrontOBJ instance
Handles only very rudimentary reading and contains no error handling!
Does not handle:
- relative indexing
- subobjects or groups
- lines, splines, beziers, etc.
"""
# parses a vertex record as either vid, vid/tid, vid//nid or vid/tid/nid
# and returns a 3-tuple where unparsed values are replaced with -1
def parse_vertex( vstr ):
vals = vstr.split('/')
vid = int(vals[0])-1
tid = int(vals[1])-1 if len(vals) > 1 and vals[1] else -1
nid = int(vals[2])-1 if len(vals) > 2 else -1
return (vid,tid,nid)
with open( filename, 'r' ) as objf:
obj = WavefrontOBJ(default_mtl=default_mtl)
obj.path = filename
cur_mat = obj.mtls.index(default_mtl)
for line in objf:
toks = line.split()
if not toks:
continue
if toks[0] == 'v':
obj.vertices.append( [ float(v) for v in toks[1:]] )
elif toks[0] == 'vn':
obj.normals.append( [ float(v) for v in toks[1:]] )
elif toks[0] == 'vt':
obj.texcoords.append( [ float(v) for v in toks[1:]] )
elif toks[0] == 'f':
poly = [ parse_vertex(vstr) for vstr in toks[1:] ]
if triangulate:
for i in range(2,len(poly)):
obj.mtlid.append( cur_mat )
obj.polygons.append( (poly[0], poly[i-1], poly[i] ) )
else:
obj.mtlid.append(cur_mat)
obj.polygons.append( poly )
elif toks[0] == 'mtllib':
obj.mtllibs.append( toks[1] )
elif toks[0] == 'usemtl':
if toks[1] not in obj.mtls:
obj.mtls.append(toks[1])
cur_mat = obj.mtls.index( toks[1] )
return obj
def save_obj( obj: WavefrontOBJ, filename: str ):
"""Saves a WavefrontOBJ object to a file
Warning: Contains no error checking!
"""
with open( filename, 'w' ) as ofile:
for mlib in obj.mtllibs:
ofile.write('mtllib {}\n'.format(mlib))
for vtx in obj.vertices:
ofile.write('v '+' '.join(['{}'.format(v) for v in vtx])+'\n')
for tex in obj.texcoords:
ofile.write('vt '+' '.join(['{}'.format(vt) for vt in tex])+'\n')
for nrm in obj.normals:
ofile.write('vn '+' '.join(['{}'.format(vn) for vn in nrm])+'\n')
if not obj.mtlid:
obj.mtlid = [-1] * len(obj.polygons)
poly_idx = np.argsort( np.array( obj.mtlid ) )
cur_mat = -1
for pid in poly_idx:
if obj.mtlid[pid] != cur_mat:
cur_mat = obj.mtlid[pid]
ofile.write('usemtl {}\n'.format(obj.mtls[cur_mat]))
pstr = 'f '
for v in obj.polygons[pid]:
# UGLY!
vstr = '{}/{}/{} '.format(v[0]+1,v[1]+1 if v[1] >= 0 else 'X', v[2]+1 if v[2] >= 0 else 'X' )
vstr = vstr.replace('/X/','//').replace('/X ', ' ')
pstr += vstr
ofile.write( pstr+'\n')
```
Ensuite pour charger un fichier OBJ, puis créer une surface, rien de plus simple:
```python
from wavefront import *
obj = load_obj( 'spot.obj')
ps.init()
ps_mesh = ps.register_surface_mesh("spot", obj.only_coordinates(), obj.only_faces() )
ps.show()
```
## Extraire les bords dans un maillage
Le code ci-dessous, à ajouter à la classe `WavefrontOBJ` vous permet d'extraire les bords dans un maillage (i.e. les arêtes incidentes à une seule face). La méthode ci-dessous retourne les sommets sur le bord, sous forme d'une liste `[]`.
```python
### return boundary vertices
def boundary_vertices( self ):
bdry_v= set()
darts = {}
faces = self.only_faces()
for f in faces:
for i in range( len(f) ):
if ( f[i] in darts ):
darts[ f[ i ] ].append( f[ (i+1)%len(f) ] )
else:
darts[ f[ i ] ] = [ f[ (i+1)%len(f) ] ]
for i, i_list in darts.items():
for j in i_list:
j_list = darts[ j ]
if ( not i in j_list ):
bdry_v.add( i )
bdry_v.add( j )
return bdry_v
```
La méthode ci-dessous retourne une liste d'arcs (arêtes orientées) sur le bord du maillage.
```python
### return boundary edges
def boundary_edges( self ):
bdry_e= set()
darts = {}
faces = self.only_faces()
for f in faces:
for i in range( len(f) ):
if ( f[i] in darts ):
darts[ f[ i ] ].append( f[ (i+1)%len(f) ] )
else:
darts[ f[ i ] ] = [ f[ (i+1)%len(f) ] ]
for i, i_list in darts.items():
for j in i_list:
j_list = darts[ j ]
if ( not i in j_list ):
bdry_e.add( (j,i) )
return bdry_e
```
La méthode ci-dessous fait la même chose mais retourne cette liste d'arcs sous la forme d'un tableau numpy à deux dimensions. Cela permet d'afficher directement ces arcs sous polyscope.
```python
### return boundary edges
def numpy_boundary_edges( self ):
edges = self.boundary_edges()
np_edges = np.zeros( ( len(edges), 2 ), np.int32 )
k = 0;
for e in edges:
np_edges[k][0]=e[0]
np_edges[k][1]=e[1]
k = k+1
return np_edges
```
Enfin, la méthode ci-dessous reconstruit la liste ordonnée des sommets sur le bord (si il y a un seul bord).
```python
def ordered_boundary( self ):
edges = self.boundary_edges()
if len(edges)==0:
return []
d = {}
first = -1
for e in edges:
d[ e[0] ] = e[1]
if first == -1:
first = e[ 0 ]
list_e = [first]
cur = d[ first ]
while cur != first:
list_e.append( cur )
cur = d[ cur ]
return list_e
```
Voilà ci-dessous un petit exemple d'utilisation.
```python=
import polyscope as ps
import numpy as np
from wavefront import *
name='bunnyhead'
obj_name=name+'.obj'
obj = load_obj( obj_name )
ps.init()
ps_mesh = ps.register_surface_mesh(name, obj.only_coordinates(), obj.only_faces() )
bdry = obj.numpy_boundary_edges()
ps_net= ps.register_curve_network("boundary", obj.only_coordinates(), bdry )
ps.show()
```
Vous pouvez récupérer la liste ordonnée des sommets sur le bord via:
```python
>>> L=obj.ordered_boundary()
>>> L
[566, 567, 547, 543, 540, 533, 529, 526, 519, 517, 515, 40, 816, 1412, 1411, 1433, 1442, 1449, 1379, 639, 57, 56, 278, 1480, 1465, 392, 119, 393, 506, 793, 227, 96, 771, 747, 309, 645, 1197, 67, 632, 613, 571, 537, 513, 483, 438, 427, 397, 386, 365, 349, 323, 1148, 1100, 1137, 1091, 1015, 1022, 1031, 1058, 1065, 1070, 1074, 1084, 322, 911, 656, 658, 1382, 249, 424, 83, 916, 250, 254, 1306, 685, 700, 718, 742, 743, 770, 792, 800, 229, 426, 1199, 436, 505, 1207, 1502, 1179, 457, 458, 1174, 446, 401, 511, 569]
```