#
##
## This file is part of pyFormex 2.0 (Mon Sep 14 12:29:05 CEST 2020)
## pyFormex is a tool for generating, manipulating and transforming 3D
## geometrical models by sequences of mathematical operations.
## Home page: http://pyformex.org
## Project page: http://savannah.nongnu.org/projects/pyformex/
## Copyright 2004-2020 (C) Benedict Verhegghe (benedict.verhegghe@ugent.be)
## Distributed under the GNU General Public License version 3 or later.
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see http://www.gnu.org/licenses/.
##
"""Operations on triangulated surfaces.
A triangulated surface is a surface consisting solely of triangles.
Any surface in space, no matter how complex, can be approximated with
a triangulated surface.
"""
import numpy as np
from pyformex import Path
from pyformex import arraytools as at
from pyformex import fileread
from pyformex import filewrite
from pyformex import geomtools
from pyformex import inertia
from pyformex import utils
from pyformex.coords import Coords
from pyformex.connectivity import Connectivity
from pyformex.mesh import Mesh
from pyformex.formex import Formex
from pyformex.geometry import Geometry
__all__ = ['TriSurface', 'fillBorder']
#
# gts commands used:
# in Debian package: stl2gts gts2stl gtscheck
# not in Debian package: gtssplit gtscoarsen gtsrefine gtssmooth gtsinside
#
# TODO: check if we can replace this with frontwalks
def _adjacencyArrays(elems, nsteps=1):
"""Create adjacency arrays for 2-node elements.
elems is a (nr,2) shaped integer array.
The result is a list of adjacency arrays, where row i of adjacency array j
holds a sorted list of the nodes that are connected to node i via a shortest
path of j elements, padded with -1 values to create an equal list length
for all nodes.
This is: [adj0, adj1, ..., adjj, ... , adjn] with n=nsteps.
Examples
>>> for a in _adjacencyArrays([[0,1],[1,2],[2,3],[3,4],[4,0]],3):
... print(a)
[[0]
[1]
[2]
[3]
[4]]
[[1 4]
[0 2]
[1 3]
[2 4]
[0 3]]
[[2 3]
[3 4]
[0 4]
[0 1]
[1 2]]
[]
"""
elems = Connectivity(elems)
if len(elems.shape) != 2 or elems.shape[1] != 2:
raise ValueError("""Expected a set of 2-node elements.""")
if nsteps < 1:
raise ValueError("""The shortest path should be at least 1.""")
# Construct table of nodes connected to each node
adj1 = elems.adjacency('n')
m = adj1.shape[0]
adj = [np.arange(m).reshape(-1, 1), adj1]
nodes = adj1
step = 2
while step <= nsteps and nodes.size > 0:
# Determine adjacent nodes
t = nodes < 0
nodes = adj1[nodes]
nodes[t] = -1
nodes = nodes.reshape((m, -1))
nodes = nodes.normalize()
# Remove nodes of lower adjacency
ladj = np.concatenate(adj[-2:], -1)
t = [np.in1d(n, l, assume_unique=True) for n, l in zip(nodes, ladj)]
t = np.asarray(t)
nodes[t] = -1
nodes = nodes.sortRows()
# Store current nodes
adj.append(nodes)
step += 1
return adj
############################################################################
# The TriSurface class
[docs]class TriSurface(Mesh):
"""A class representing a triangulated 3D surface.
A triangulated surface is a surface consisting of a collection of
triangles. The TriSurface is subclassed from :class:`Mesh` with a
fixed plexitude 3.
The surface contains `ntri` triangles and `nedg` edges. Each triangle
has 3 vertices with 3 coordinates. The total number of vertices is
`ncoords`.
The TriSurface can be initialized from one of the following sets of data:
- an (ntri,3,3) shaped array of floats
- a Formex with plexitude 3
- a Mesh with plexitude 3
- an (ncoords,3) float array of vertex coordinates and
an (ntri,3) integer array of vertex numbers
- an (ncoords,3) float array of vertex coordinates,
an (nedg,2) integer array of vertex numbers,
an (ntri,3) integer array of edges numbers.
Additionally, a keyword argument ``prop`` may be provided to
specify property values, as in :class:`Mesh`.
"""
_exclude_members_ = ['intersectionWithLines']
def __init__(self, *args, **kargs):
"""Create a new surface."""
self._areas = self._fnormals = None
self.adj = None
if hasattr(self, 'edglen'):
del self.edglen
if len(args) == 0:
Mesh.__init__(self, [], [], None, 'tri3')
return # an empty surface
if len(args) == 1:
# argument should be a suitably structured geometry object
# TriSurface, Mesh, Formex, Coords, ndarray, ...
a = args[0]
if isinstance(a, Mesh):
if a.nplex() != 3 or a.elName() != 'tri3':
raise ValueError(
"Only meshes with plexitude 3 and eltype 'tri3' "
"can be converted to TriSurface!")
Mesh.__init__(self, a.coords, a.elems, a.prop, 'tri3')
else:
if not isinstance(a, Formex):
# something that can be converted to a Formex
try:
a = Formex(a)
except ValueError:
raise ValueError(
"Can not convert objects of type "
f"{type(a)} to TriSurface!")
# now a is a Formex
if a.nplex() != 3:
raise ValueError("Expected an object with plexitude 3!")
coords, elems = a.coords.fuse()
Mesh.__init__(self, coords, elems, a.prop, 'tri3')
else:
# arguments are (coords,elems) or (coords,edges,faces)
coords = Coords(args[0])
if len(coords.shape) != 2:
raise ValueError("Expected a 2-dim coordinates array")
if len(args) == 2:
# arguments are (coords,elems)
elems = Connectivity(args[1], nplex=3)
Mesh.__init__(self, coords, elems, None, 'tri3')
elif len(args) == 3:
# arguments are (coords,edges,faces)
edges = Connectivity(args[1], nplex=2)
if edges.size > 0 and edges.max() >= coords.shape[0]:
raise ValueError("Some vertex number is too high")
faces = Connectivity(args[2], nplex=3)
if faces.max() >= edges.shape[0]:
raise ValueError("Some edge number is too high")
elems = faces.combine(edges)
Mesh.__init__(self, coords, elems, None, 'tri3')
# since we have the extra data available, keep them
self.edges = edges
self.elem_edges = faces
else:
raise RuntimeError("Too many positional arguments")
if 'prop' in kargs:
self.setProp(kargs['prop'])
def __setstate__(self, state):
"""Set the object from serialized state.
This allows to read back old pyFormex Project files where the
Surface class did not set an element type.
"""
if 'areas' in state:
state['_areas'] = state['areas']
del state['areas']
self.__dict__.update(state)
self.eltype = 'tri3'
#######################################################################
#
# Return information about a TriSurface
#
[docs] def nedges(self):
"""Return the number of edges of the TriSurface."""
return self.getEdges().shape[0]
[docs] def nfaces(self):
"""Return the number of faces of the TriSurface."""
return self.getElemEdges().shape[0]
[docs] def vertices(self):
"""Return the coordinates of the nodes of the TriSurface."""
return self.coords
[docs] def shape(self):
"""Return the number of points, edges, faces of the TriSurface."""
return self.ncoords(), self.nedges(), self.nfaces()
[docs] def getElemEdges(self):
"""Get the faces' edge numbers."""
if self.elem_edges is None:
self.elem_edges, self.edges = self.elems.insertLevel(1)
return self.elem_edges
#######################################################################
#
# Operations that change the TriSurface itself
#
# Make sure that you know what you're doing if you use these
#
#
# Changes to the geometry should by preference be done through the
# __init__ function, to ensure consistency of the data.
# Convenience functions are defined to change some of the data.
#
[docs] def setCoords(self, coords):
"""Change the coords."""
self.__init__(coords, self.elems, prop=self.prop)
return self
[docs] def setElems(self, elems):
"""Change the elems."""
self.__init__(self.coords, elems, prop=self.prop)
[docs] def setEdgesAndFaces(self, edges, faces):
"""Change the edges and faces."""
self.__init__(self.coords, edges, faces, prop=self.prop)
[docs] def append(self, S):
"""Merge another surface with self.
This just merges the data sets, and does not check
whether the surfaces intersect or are connected!
This is intended mostly for use inside higher level functions.
"""
coords = np.concatenate([self.coords, S.coords])
elems = np.concatenate([self.elems, S.elems+self.ncoords()])
## What to do if one of the surfaces has properties, the other one not?
## The current policy is to use zero property values for the Surface
## without props
prop = None
if self.prop is not None or S.prop is not None:
if self.prop is None:
self.prop = np.zeros(shape=self.nelems(), dtype=at.Int)
if S.prop is None:
p = np.zeros(shape=S.nelems(), dtype=at.Int)
else:
p = S.prop
prop = np.concatenate((self.prop, p))
self.__init__(coords, elems, prop=prop)
#######################################################################
#
# read and write
#
[docs] @classmethod
def read(clas, fn, ftype=None):
"""Read a surface from file.
If no file type is specified, it is derived from the filename
extension.
Currently supported file types:
- .off
- .gts
- .stl (ASCII or BINARY)
- .neu (Gambit Neutral)
- .smesh (Tetgen)
Compressed (gzip or bzip2) files are also supported. Their names
should be the normal filename with '.gz' or '.bz2' appended.
These files are uncompressed on the fly during the reading and
the uncompressed versions are deleted after reading.
The file type can be specified explicitely to handle file names
where the extension does not directly specify the file type.
"""
fn = Path(fn)
if ftype is None:
ftype, compr = fn.ftype_compr()
else:
ftype, compr = Path('a.'+ftype).ftype_compr()
if ftype == 'off':
data = fileread.readOFF(fn)
elif ftype == 'gts':
data = fileread.readGTS(fn)
elif ftype == 'stl':
coords, normals, color = fileread.readSTL(fn)
S = TriSurface(coords)
# S.addField('elemc', normals, '_fnormals_')
# S.setNormals(S.getField('_fnormals_').convert('elemn'))
if color:
S.attrib(color=color)
return S
# The remainder should probably disappear
# The following do not support compression yet
elif ftype == 'smesh' and not compr:
from pyformex.plugins import tetgen
data = tetgen.readSurface(fn)
elif ftype == 'neu' and not compr:
data = fileread.read_gambit_neutral(fn)
elif ftype in ['vtp', 'vtk'] and not compr:
return read_vtk_surface(fn)
else:
raise f"Unknown TriSurface type, cannot read file {fn}"
return TriSurface(*data)
[docs] def write(self, fname, ftype=None, color=None):
"""Write the surface to file.
If no filetype is given, it is deduced from the filename extension.
If the filename has no extension, the 'off' file type is used.
For a file with extension 'stl', the ftype may be 'stla' or 'stlb'
to force ascii or binary STL format.
The color is only useful for 'stlb' format.
"""
fname = Path(fname)
if ftype is None:
ftype, compr = fname.ftype_compr()
else:
ftype, compr = Path('a.'+ftype).ftype_compr()
if compr and ftype in ['smesh', 'vtp', 'vtk']:
raise ValueError("Compressed surface export is not active (yet)")
print(f"Writing surface to file {fname} ({ftype})")
if ftype == 'pgf':
Geometry.write(self, fname)
elif ftype == 'gts':
filewrite.writeGTS(fname, self)
print("Wrote {} vertices, {} edges, {} faces".format(*self.shape()))
elif ftype in ['stl', 'stla', 'stlb', 'off', 'smesh', 'vtp', 'vtk']:
if ftype in ['stl', 'stla']:
filewrite.writeSTL(fname, self.coords[self.elems], binary=False)
elif ftype == 'stlb':
filewrite.writeSTL(fname, self.coords[self.elems], binary=True,
color=color)
elif ftype == 'off':
filewrite.writeOFF(fname, self)
elif ftype == 'smesh':
from pyformex.plugins import tetgen
tetgen.writeSurface(fname, self.coords, self.elems)
elif ftype == 'vtp' or ftype == 'vtk':
from pyformex.plugins import vtk_itf
vtk_itf.writeVTP(fname, self, checkMesh=False)
print(f"Wrote {self.ncoords()} vertices, {self.nelems()} elems")
else:
print(f"Cannot save TriSurface as file {fname}")
####################### TriSurface Data ######################
[docs] def avgVertexNormals(self):
"""Compute the average normals at the vertices."""
return geomtools.averageNormals(self.coords, self.elems, atNodes=True)
[docs] def areaNormals(self):
"""Compute the area and normal vectors of the surface triangles.
The normal vectors are normalized.
The area is always positive.
The values are returned and saved in the object.
"""
if self._areas is None or self._fnormals is None:
self._areas, self._fnormals = geomtools.areaNormals(self.coords[self.elems])
return self._areas, self._fnormals
[docs] def areas(self):
"""Return the areas of all facets"""
return self.areaNormals()[0]
[docs] def volume(self):
"""Return the enclosed volume of the surface.
This will only be correct if the surface is a closed manifold.
"""
x = self.coords[self.elems]
return inertia.surface_volume(x).sum()
[docs] def volumeInertia(self, density=1.0):
"""Return the inertia properties of the enclosed volume of the surface.
The surface should be a closed manifold and is supposed to be
the border of a volume of constant density 1.
Returns an :class:`inertia.Inertia` instance with attributes
- `mass`: the total mass (float)
- `ctr`: the center of mass: float (3,)
- `tensor`: the inertia tensor in the central axes: shape (3,3)
This will only be correct if the surface is a closed manifold.
See :meth:`inertia` for the inertia of the surface.
Example:
>>> from pyformex.simple import sphere
>>> S = sphere(8)
>>> I = S.volumeInertia()
>>> print(I.mass) # doctest: +ELLIPSIS
4.1526...
>>> print(I.ctr)
[ 0. 0. 0.]
>>> print(I.tensor)
[[ 1.65 0. -0. ]
[ 0. 1.65 -0. ]
[-0. -0. 1.65]]
"""
x = self.coords[self.elems]
V, C, I = inertia.surface_volume_inertia(x) # noqa: E741
I = inertia.Tensor(I) # noqa: E741
I = inertia.Inertia(I, mass=V, ctr=C) # noqa: E741
I.mass *= density
I *= density # noqa: E741
return I
[docs] def inertia(self, volume=False, density=1.0):
"""Return inertia related quantities of the surface.
This computes the inertia properties of the centroids of the
triangles, using the triangle area as a weight. The result is
therefore different from self.coords.inertia() and usually better
suited for the surface, especially if the triangle areas differ a lot.
Returns a tuple with the center of gravity, the principal axes of
inertia, the principal moments of inertia and the inertia tensor.
See also :meth:`volumeInertia`.
"""
if volume:
return self.volumeInertia(density=density)
else:
I = self.centroids().inertia(mass=self.areas()) # noqa: E741
I.mass *= density
I *= density # noqa: E741
return I
[docs] def curvature(self, neighbours=1):
"""Calculate curvature parameters at the nodes.
Algorithms based on Dong and Wang 2005; Koenderink and Van Doorn 1992.
This uses the nodes that are connected to the node via a shortest
path of 'neighbours' edges.
Eight values are returned: the Gaussian and mean curvature, the
shape index, the curvedness, the principal curvatures and the
principal directions.
"""
# calculate n-ring neighbourhood of the nodes (n=neighbours)
adj = _adjacencyArrays(self.getEdges(), nsteps=neighbours)[-1]
adjNotOk = adj<0
# remove the nodes that have less than three adjacent nodes,
adjNotOk[(adj>=0).sum(-1) <= 2] = True
# calculate unit length average normals at the nodes p
# a weight 1/|gi-p| could be used (gi=center of the face fi)
p = self.coords
n = geomtools.averageNormals(self.coords, self.elems, atNodes=True)
# double-precision: this will allow us to check the sign of the angles
# BV: why is double precision needed to check a sign???
p = p.astype(np.float64)
n = n.astype(np.float64)
vp = p[adj] - p[:, np.newaxis]
vn = n[adj] - n[:, np.newaxis]
# where adjNotOk, set vectors = [0.,0.,0.]
# this will result in NaN values
vp[adjNotOk] = 0.
vn[adjNotOk] = 0.
# calculate unit length projection of vp onto the tangent plane
t = geomtools.projectionVOP(vp, n[:, np.newaxis])
t = at.normalize(t)
# calculate normal curvature
with np.errstate(divide='ignore', invalid='ignore'):
k = at.dotpr(vp, vn) / at.dotpr(vp, vp)
# calculate maximum normal curvature and corresponding coordinate system
try:
imax = np.nanargmax(k, -1)
kmax = k[np.arange(len(k)), imax]
tmax = t[np.arange(len(k)), imax]
except Exception:
# bug with numpy.nanargmax: cannot convert float NaN to integer
kmax = np.resize(np.NaN, (k.shape[0]))
tmax = np.resize(np.NaN, (t.shape[0], 3))
w = ~(np.isnan(k).all(1))
imax = np.nanargmax(k[w], -1)
kmax[w] = k[w, imax]
tmax[w] = t[w, imax]
tmax1 = tmax
tmax2 = np.cross(n, tmax1)
tmax2 = at.normalize(tmax2)
# calculate angles (tmax1,t)
theta, rot = geomtools.rotationAngle(
np.repeat(tmax1[:, np.newaxis], t.shape[1], 1),
t, angle_spec=at.RAD)
theta = theta.reshape(t.shape[:2])
rot = rot.reshape(t.shape)
# check the sign of the angles
d = at.dotpr(rot, n[:, np.newaxis])/(
# divide by length for round-off errors
at.length(rot)*at.length(n)[:, np.newaxis])
cw = np.isclose(d, [-1.])
theta[cw] = -theta[cw]
# calculate coefficients
a = kmax
a11 = np.nansum(np.cos(theta)**2*np.sin(theta)**2, -1)
a12 = np.nansum(np.cos(theta)*np.sin(theta)**3, -1)
a22 = np.nansum(np.sin(theta)**4, -1)
a13 = np.nansum((k-a[:, np.newaxis] * np.cos(theta)**2)
* np.cos(theta)*np.sin(theta), -1)
a23 = np.nansum((k-a[:, np.newaxis]*np.cos(theta)**2)*np.sin(theta)**2, -1)
denom = (a11*a22-a12**2)
b = (a13*a22-a23*a12)/denom
c = (a11*a23-a12*a13)/denom
# calculate the Gaussian and mean curvature
K = a*c-b**2/4
H = (a+c)/2
# calculate the principal curvatures and principal directions
k1 = H+np.sqrt(H**2-K)
k2 = H-np.sqrt(H**2-K)
theta0 = 0.5*np.arcsin(b/(k2-k1))
w = np.apply_along_axis(np.isclose, 0, -b, 2*(k2-k1)*np.cos(theta0)*np.sin(theta0))
theta0[w] = np.pi-theta0[w]
e1 = np.cos(theta0)[:, np.newaxis]*tmax1 +np.sin(theta0)[:, np.newaxis]*tmax2
e2 = np.cos(theta0)[:, np.newaxis]*tmax2 -np.sin(theta0)[:, np.newaxis]*tmax1
# calculate the shape index and curvedness
S = 2./np.pi * np.arctan((k1+k2)/(k1-k2))
C = np.square((k1**2 + k2**2) / 2)
return K, H, S, C, k1, k2, e1, e2
return curv
[docs] def surfaceType(self):
"""Check whether the TriSurface is a manifold and if it's closed."""
ncon = self.nEdgeConnected()
maxcon = ncon.max()
mincon = ncon.min()
manifold = maxcon == 2
closed = mincon == 2
return manifold, closed, mincon, maxcon
[docs] def borderEdges(self):
"""Detect the border elements of TriSurface.
The border elements are the edges having less than 2 connected elements.
Returns True where edge is on the border.
"""
return self.nEdgeConnected() <= 1
[docs] def borderEdgeNrs(self):
"""Returns the numbers of the border edges."""
return np.where(self.nEdgeConnected() <= 1)[0]
[docs] def borderNodeNrs(self):
"""Detect the border nodes of TriSurface.
The border nodes are the vertices belonging to the border edges.
Returns a list of vertex numbers.
"""
border = self.getEdges()[self.borderEdgeNrs()]
return np.unique(border)
####### MANIFOLD #################
[docs] def isManifold(self):
"""Check whether the TriSurface is a manifold.
A surface is a manifold if a small sphere exists that cuts the surface
to a surface that can continously be deformed to an open disk.
"""
return self.surfaceType()[0]
# TODO: We should probably add optional sorting (# connections) here
# TODO: this could be moved into Mesh.nonManifoldEdges if it works
# for all level 2 Meshes
[docs] def nonManifoldEdges(self):
"""Return the non-manifold edges.
Non-manifold edges are edges having more than two triangles
connected to them.
Returns the indices of the non-manifold edges in a TriSurface.
"""
return np.where(self.nEdgeConnected()>2)[0]
[docs] def nonManifoldEdgesFaces(self):
"""Return the non-manifold edges and faces.
Non-manifold edges are edges that are connected to more than
two faces.
Returns
-------
edges: list of int
The list of non-manifold edges.
faces: list of int
The list of faces connected to any of the non-manifold edges.
"""
conn = self.edgeConnections()
ed = (conn!=-1).sum(axis=1)>2
fa = np.unique(conn[ed])
return np.arange(len(ed))[ed], fa[fa!=-1]
[docs] def isClosedManifold(self):
"""Check whether the TriSurface is a closed manifold."""
stype = self.surfaceType()
return stype[0] and stype[1]
[docs] def isConvexManifold(self):
"""Check whether the TriSurface is a convex manifold."""
return self.isManifold() and self.edgeSignedAngles().min()>=0.
[docs] def removeNonManifold(self):
"""Remove the non-manifold edges.
Removes the non-manifold edges by iteratively applying
:meth:`removeDuplicate` and :meth:`collapseEdge` until no edge
has more than two connected triangles.
Returns the reduced surface.
"""
S = self.removeDuplicate()
non_manifold_edges = self.nonManifoldEdges()
while non_manifold_edges.any():
print(f"# nonmanifold edges: {len(non_manifold_edges)}")
maxcon = S.nEdgeConnected().max()
wmax = np.where(S.nEdgeConnected()==maxcon)[0]
S = S.collapseEdge(wmax[0])
S = S.removeDuplicate()
non_manifold_edges = S.nonManifoldEdges()
return S
###### BORDER ######################
[docs] def checkBorder(self):
"""Return the border of TriSurface.
Returns a list of connectivity tables. Each table holds the
subsequent line segments of one continuous contour of the border
of the surface.
"""
border = self.getEdges()[self.borderEdges()]
if len(border) > 0:
return border.chained()
else:
return []
[docs] def border(self, compact=True):
"""Return the border(s) of TriSurface.
The complete border of the surface is returned as a list
of plex-2 Meshes. Each Mesh constitutes a continuous part
of the border. By default, the Meshes are compacted.
Setting compact=False will return all Meshes with the full
surface coordinate sets. This is usefull for filling the
border and adding to the surface.
"""
ML = [Mesh(self.coords, e) for e in self.checkBorder()]
if compact:
ML = [M.compact() for M in ML]
return ML
[docs] def fillBorder(self, method='radial', dir=None, compact=True):
"""Fill the border areas of a surface to make it closed.
Returns a list of surfaces, each of which fills a singly connected
part of the border of the input surface. Adding these surfaces to
the original will create a closed surface. The surfaces will have
property values set above those used in the parent surface.
If the surface is already closed, an empty list is returned.
There are three methods: 'radial', 'planar' and 'border',
corresponding to the methods of the :func:`fillBorder` function.
"""
if self.prop is None:
mprop = 1
else:
mprop = self.prop.max()+1
return [fillBorder(b, method, dir).setProp(mprop+i)
for i, b in enumerate(self.border(compact=compact))]
[docs] def close(self, method='radial', dir=None):
"""This method needs documentation!!!!"""
border = self.fillBorder(method, dir, compact=method=='radial')
if method == 'radial':
return self.concatenate([self]+border)
else:
elems = np.concatenate([m.elems for m in [self]+border], axis=0)
if self.prop is None:
prop = np.zeros(shape=self.nelems(), dtype=at.Int)
else:
prop = self.prop
prop = np.concatenate([prop] + [m.prop for m in border])
return TriSurface(self.coords, elems, prop=prop)
[docs] def edgeCosAngles(self, return_mask=False):
"""Return the cos of the angles over all edges.
The surface should be a manifold (max. 2 elements per edge).
Edges adjacent to only one element get cosangles = 1.0.
If return_mask == True, a second return value is a boolean array
with the edges that connect two faces.
As a side effect, this method also sets the area, normals,
elem_edges and edges attributes.
"""
# get connections of edges to faces
conn = self.getElemEdges().inverse(expand=True)
## # BV: The following gives the same results, but does not
## # guarantee the edges attribute to be set
## conn1 = self.edgeConnections()
## diff = (conn1 != conn).sum()
## if diff > 0:
## print "edgeConnections",conn1
## print "getElemEdges.inverse",conn
# Bail out if some edge has more than two connected faces
if conn.shape[1] > 2:
raise RuntimeError("The TriSurface is not a manifold")
# get normals on all faces
n = self.areaNormals()[1]
# Flag edges that connect two faces
conn2 = (conn >= 0).sum(axis=-1) == 2
# get adjacent facet normals for 2-connected edges
n = n[conn[conn2]]
# compute cosinus of angles over 2-connected edges
cosa = at.dotpr(n[:, 0], n[:, 1])
# Initialize cosangles to all 1. values
cosangles = np.ones((conn.shape[0],))
# Fill in the values for the 2-connected edges
cosangles[conn2] = cosa
# Clip to the -1...+1. range
cosangles = cosangles.clip(min=-1., max=1.)
# Return results
if return_mask:
return cosangles, conn2
else:
return cosangles
[docs] def edgeAngles(self):
"""Return the angles over all edges (in degrees).
It is the angle (0 to 180) between 2 face normals.
"""
return at.arccosd(self.edgeCosAngles())
[docs] def edgeSignedAngles(self, return_mask=False):
"""Return the signed angles over all edges (in degrees).
It is the angle (-180 to 180) between 2 face normals.
Positive/negative angles are associated to convexity/concavity
at that edge.
The border edges attached to one triangle have angle 0.
NB: The sign of the angle is relevant if the surface has fixed normals.
Should this check be done?
"""
# get connections of edges to faces
conn = self.getElemEdges().inverse(expand=True)
if conn.shape[1] > 2:
raise RuntimeError("The TriSurface is not a manifold")
# get normals on all faces
n = self.areaNormals()[1]
# Flag edges that connect two faces
conn2 = (conn >= 0).sum(axis=-1) == 2
# get adjacent facet normals for 2-connected edges
n = n[conn[conn2]]
edg = self.coords[self.getEdges()]
edg = edg[:, 1]-edg[:, 0]
ang = geomtools.rotationAngle(n[:, 0], n[:, 1], m=edg[conn2],
angle_spec=at.DEG)
# Initialize signed angles to all 0. values
sangles = np.ones((conn.shape[0],))
sangles[conn2] = ang
# Clip to the -180...+180. range
sangles = sangles.clip(min=-180., max=180.)
# Return results
if return_mask:
return sangles, conn2
else:
return sangles
[docs] def edgeLengths(self):
"""Returns the lengths of all edges
Returns an array with the length of all the edges in the
surface. As a side effect, this stores the connectivities
of the edges to nodes and the elements to edges in the
attributes edges, resp. elem_edges.
"""
edg = self.coords[self.getEdges()]
return at.length(edg[:, 1]-edg[:, 0])
def _compute_data(self):
"""Compute data for all edges and faces."""
if hasattr(self, 'edglen'):
return
self.areaNormals()
self.edglen = self.edgeLengths()
self.facedg = self.edglen[self.getElemEdges()]
self.peri = self.facedg.sum(axis=-1)
self.edgmin = self.facedg.min(axis=-1)
self.edgmax = self.facedg.max(axis=-1)
self.altmin = 2*self._areas / self.edgmax
self.aspect = self.edgmax/self.altmin
qual_equi = np.sqrt(np.sqrt(3.)) / 6.
self.qual = np.sqrt(self._areas) / self.peri / qual_equi
[docs] def perimeters(self):
"""Compute the perimeters of all triangles."""
self._compute_data()
return self.peri
[docs] def quality(self):
"""Compute a quality measure for the triangle schapes.
The quality of a triangle is defined as the ratio of the square
root of its surface area to its perimeter relative to this same
ratio for an equilateral triangle with the same area. The quality
is then one for an equilateral triangle and tends to zero for a
very stretched triangle.
"""
self._compute_data()
return self.qual
[docs] def aspectRatio(self):
"""Return the apect ratio of the triangles of the surface.
The aspect ratio of a triangle is the ratio of the longest edge
over the smallest altitude of the triangle.
Equilateral triangles have the smallest edge ratio (2 over square root 3).
"""
self._compute_data()
return self.aspect
[docs] def smallestAltitude(self):
"""Return the smallest altitude of the triangles of the surface."""
self._compute_data()
return self.altmin
[docs] def longestEdge(self):
"""Return the longest edge of the triangles of the surface."""
self._compute_data()
return self.edgmax
[docs] def shortestEdge(self):
"""Return the shortest edge of the triangles of the surface."""
self._compute_data()
return self.edgmin
[docs] def stats(self):
"""Return a text with full statistics."""
bbox = self.bbox()
manifold, closed, mincon, maxcon = self.surfaceType()
self._compute_data()
area = self.area()
qual = self.quality()
s = f"""
Size: {self.ncoords()} vertices, {self.nedges()} edges and {self.nfaces()} faces
Bounding box: min {bbox[0]}, max {bbox[1]}
Minimal/maximal number of connected faces per edge: {mincon}/{maxcon}
Surface is{'' if manifold else ' not'} a{' closed' if closed else ''} manifold
Facet area: min { self._areas.min()}; mean {self._areas.mean()}; max {self._areas.max()}
Edge length: min {self.edglen.min()}; mean {self.edglen.mean()}; max {self.edglen.max()}
Shortest altitude: {self.altmin.min()}; largest aspect ratio: {self.aspect.max()}
Quality: {qual.min()} .. {qual.max()}
"""
if manifold:
angles = self.edgeAngles()
# getAngles is currently removed
# vangles = self.getAngles()
if closed:
volume = self.volume()
s += (f"Angle between adjacent facets: min: {angles.min()}; "
f"mean: {angles.mean()}; max: {angles.max()}\n")
s += f"Total area: {area}; "
if manifold and closed:
s += f"Enclosed volume: {volume}"
else:
s += "No volume (not a closed manifold!)"
return s
[docs] def distanceOfPoints(self, X, return_points=False):
"""Find the distances of points X to the TriSurface.
The distance of a point is either:
- the closest perpendicular distance to the facets;
- the closest perpendicular distance to the edges;
- the closest distance to the vertices.
X is a (nX,3) shaped array of points.
If return_points = True, a second value is returned: an array with
the closest (foot)points matching X.
"""
from pyformex.timer import Timer
t = Timer()
# distance from vertices
ind, dist = geomtools.closest(X, self.coords, return_dist=True)
if return_points:
points = self.coords[ind]
print(f"Vertex distance: {t.seconds(True)} seconds")
# distance from edges
Ep = self.coords[self.getEdges()]
res = geomtools.edgeDistance(X, Ep, return_points) # OKpid, OKdist, (OKpoints)
okE, distE = res[:2]
closer = distE < dist[okE]
if closer.size > 0:
dist[okE[closer]] = distE[closer]
if return_points:
points[okE[closer]] = res[2][closer]
print(f"Edge distance: {t.seconds(True)} seconds")
# distance from faces
Fp = self.coords[self.elems]
res = geomtools.faceDistance(X, Fp, return_points) # OKpid, OKdist, (OKpoints)
okF, distF = res[:2]
closer = distF < dist[okF]
if closer.size > 0:
dist[okF[closer]] = distF[closer]
if return_points:
points[okF[closer]] = res[2][closer]
print(f"Face distance: {t.seconds(True)} seconds")
if return_points:
return dist, points
else:
return dist
[docs] def degenerate(self):
"""Return a list of the degenerate faces according to area and normals.
A face is degenerate if its surface is less or equal to zero or the
normal has a nan.
Returns the list of degenerate element numbers in a sorted array.
"""
return geomtools.degenerate(*self.areaNormals())
[docs] def removeDegenerate(self, compact=False):
"""Remove the degenerate elements from a TriSurface.
Returns a TriSurface with all degenerate elements removed.
By default, the coords set is unaltered and will still contain
all points, even ones that are no longer connected to any element.
To reduce the coordinate set, set the compact argument to True
or use the :meth:`compact` method afterwards.
"""
return self.cselect(self.degenerate(), compact=compact)
[docs] def collapseEdge(self, edg):
"""Collapse an edge in a TriSurface.
Collapsing an edge removes the triangles connected to the edge
and replaces the two vertices of the edge with a single one,
placed at the center of the edge.
Triangles connected to one of the edge vertices, will
become connected to the new vertex.
Parameters
----------
edg: int
The index of the edg to be removed. This is an index in the
array of edges as returned by :meth:`getElemEdges`.
Returns
-------
TriSurface
An almost equivalent surface with the specified edge removed.
"""
# remove the elements connected to the collapsing edge
invee = self.getElemEdges().inverse(expand=True)
els = invee[edg]
els = els[els>=0]
keep = at.complement(els, n=self.nelems())
elems = self.elems[keep]
prop = self.prop
if prop is not None:
prop = prop[keep]
# replace the first node with the mean of the two
# TODO: If one vertex is on the border and the other is not,
# use the border vertex coordinates instead of the mean
node0, node1 = self.edges[edg]
elems[elems==node1] = node0
coords = self.coords.copy()
coords[node0] = 0.5 * (coords[node0] + coords[node1])
return TriSurface(coords, elems, prop=prop).compact()
############## Transform surface #############################
# All transformations return a new surface
[docs] def offset(self, distance=1.):
"""Offset a surface with a certain distance.
All the nodes of the surface are translated over a specified distance
along their normal vector.
"""
n = self.avgVertexNormals()
coordsNew = self.coords + n*distance
return TriSurface(coordsNew, self.elems, prop=self.prop)
[docs] def dualMesh(self, method='median'):
"""Return the dual mesh of a compacted triangulated surface.
It creates a new triangular mesh where all triangles with prop `p`
represent the dual mesh region around the original surface node `p`.
For more info, see http://users.led-inc.eu/~phk/mesh-dualmesh.html.
- `method`: 'median' or 'voronoi'.
Returns:
- `method` = 'median': the Median dual mesh and the area of the region
around each node. The sum of the node-based areas is equal to the
original surface area.
- `method` = 'voronoi': the Voronoi polyeders and a None.
"""
if self.ncoords()!=self.compact().ncoords():
raise ValueError("Expected a compacted surface")
Q = self.convert('quad4', fuse=False)
if method == 'voronoi':
from pyformex.geomtools import triangleCircumCircle
Q.coords[-self.nelems():] = triangleCircumCircle(
self.coords[self.elems], bounding=False)[1]
nconn = Q.nodeConnections()[np.arange(self.ncoords())]
p = np.zeros(Q.nelems(), dtype=int)
for i, conn in enumerate(nconn):
p[conn[conn>-1]]=i
Q = Q.setProp(p)
if method == 'voronoi':
return Q, None
nodalAreas = np.asarray([Q.selectProp(i).area()
for i in range(len(Q.propSet()))])
return Q, nodalAreas
################## Partitioning a surface #############################
[docs] def featureEdges(self, angle=60.):
"""Return the feature edges of the surface.
Feature edges are edges that are prominent features of the geometry.
They are either border edges or edges where the normals on the two
adjacent triangles differ more than a given angle.
The non feature edges then represent edges on a rather smooth surface.
Parameters
----------
angle: float
The minimum value of the angle (in degrees) between the normals on
two adjacent triangles in order for the edge to be considered a
feature edge.
Returns
-------
bool array
An array with shape (nedg,) where the feature angles
are marked True.
Notes
-----
As a side effect, this also sets the `elem_edges` and `edges`
attributes, which can be used to get the edge data with the same
numbering as used in the returned mask. Thus, the following
constructs a Mesh with the feature edges of a surface S::
p = S.featureEdges()
Mesh(S.coords,S.edges[p])
"""
# Get the edge angles
cosangles, conn2 = self.edgeCosAngles(return_mask=True)
# initialize all edges as features
feature = np.ones((self.edges.shape[0],), dtype=bool)
# unmark edges with small angle
feature[conn2] = cosangles[conn2] <= at.cosd(angle)
return feature
[docs] def partitionByAngle(self, angle=60., sort='number'):
"""Partition the surface by splitting it at sharp edges.
The surface is partitioned in parts in which all elements can be
reached without ever crossing a sharp edge angle. More precisely,
any two triangles will belong to the same part if the can be connected
by a line in the surface that does not cross an edge between
two elements having their normals differ more than the specified
angle.
Parameters
----------
angle: float
The minimum value of the angle (in degrees) between the normals on
two adjacent triangles in order for the edge to be considered a
sharp edge.
sort: str
Defines how the resulting parts are sorted (by assigning them
increasing part numbers). The following sort criteria are currently
defined (any other value will return the parts unsorted):
- 'number': sort in decreasing order of the number of triangles
in the part. This is the default.
- 'area': sort according to decreasing surface area of the part.
Returns
-------
int array
An int array specifying for each triangle to which part it belongs.
Values are in the range 0..nparts.
Notes
-----
In order for the operation to be non-trivial, the specified edges,
possibly together with (parts of) the border, should form one or
more closed loops.
Beware that the existence of degenerate elements may cause
unexpected results. If unsure, use the :meth:`removeDegenerate`
method first to remove those elements.
"""
feat = self.featureEdges(angle=angle)
return self.partitionByCurve(feat, sort)
# This may replace CutWithPlane after it has been proved stable
# and has been expanded to multiple planes
[docs] def cutWithPlane1(self, p, n, side='', return_intersection=False, atol=0.):
"""Cut a surface with a plane.
Cut the surface with a plane defined by a point p and normal n.
.. warning:: This is experimental and may not work correctly.
Parameters
----------
p: float :term:`array_like` (3,)
A point in the cutting plane.
n: float :term:`array_like` (3,)
The normal vector to the cutting plane.
side: '' | '+' | '-'
Selects the returned parts. Default ('') is to return a tuple
of two surfaces, with the parts at the positive,
resp. negative side of the plane, as defined by the normal vector.
If a '+' or '-' is specified, only the corresponding part
is returned.
Returns
-------
Spos: TriSurface, optional
The part of the surfacec at the positive side of thr plane (p,n).
Only returned if side is '' or '+'.
Sneg: TriSurface, optional
The part of the surfacec at the negative side of thr plane (p,n).
Only returned if side is '' or '-'.
The returned surfaces have their normals fixed wherever possible.
Prop values are set containing the triangle number in the
original surface from which the elements resulted.
"""
def finalize(Sp, Sn, I):
# Result
res = []
if side in '+':
Sp = Sp.compact() # .fixNormals()
res.append(Sp)
if side in '-':
Sn = Sn.compact() # .fixNormals()
res.append(Sn)
if return_intersection:
res.append(I)
if len(res) == 1:
res = res[0]
else:
res = tuple(res)
return res
from pyformex.formex import _sane_side, _select_side
side = _sane_side(side)
p = at.checkArray(p, size=3, kind='f', allow='i').reshape(3)
n = at.checkArray(n, size=3, kind='f', allow='i').reshape(3)
# Make sure we inherit element number
save_prop = self.prop
self.prop = np.arange(self.elems.shape[0])
# Compute distance to plane of all vertices
d = self.distanceFromPlane(p, n)
p_pos = d > 0.
p_neg = d < 0.
p_in = ~(p_pos+p_neg)
p_posin = p_pos + p_in
p_negin = p_neg + p_in
# Reduce the surface to the part intersecting with the plane:
# Remember triangles with all vertices at same side
# Elements completely in the plane end up in both parts.
# BV: SHOULD WE CHANGE THIS???
# TODO: put them only in negative?, as the volume is at the
# negative side of the elements.
all_p = p_posin[self.elems].all(axis=-1)
all_n = p_negin[self.elems].all(axis=-1)
S = self.cclip(all_p+all_n, compact=False) # DOES THIS COMPACT? NO
Sp = self.clip(all_p, compact=False)
Sn = self.clip(all_n, compact=False)
# Restore properties
self.prop = save_prop
# If there is no intersection, we're done
# (we might have cut right along facet edges!)
if S.nelems() == 0:
res = _select_side(side, [Sp, Sn])
return res
#######################
# Cut S with the plane.
#######################
# First define facets in terms of edges
coords = S.coords
edg = S.getEdges()
fac = S.getElemEdges()
ele = S.elems
# Get the edges intersecting with the plane: 1 up and 1 down vertex
d_edg = d[edg]
edg_1_up = (d_edg > 0.).sum(axis=1) == 1
edg_1_do = (d_edg < 0.).sum(axis=1) == 1
cutedg = edg_1_up * edg_1_do
ind = np.where(cutedg)[0]
if ind.size == 0:
raise ValueError("This really should not happen!")
# Compute the intersection points
M = Mesh(S.coords, edg[cutedg])
x = geomtools.intersectionPointsSWP(
M.toFormex().coords, p, n, mode='pair', return_all=True,
atol=atol).reshape(-1, 3)
# Create inverse index to lookup the point using the edge number
rev = at.inverseUniqueIndex(ind) + coords.shape[0]
# Concatenate the coords arrays
coords = coords.concatenate([coords, x])
# For each triangle, compute the number of cutting edges
cut = cutedg[fac]
ncut = cut.sum(axis=1)
if (ncut < 1).any() or (ncut > 2).any():
# Maybe we should issue a warning and ignore these cases?
print("NCUT: ", ncut)
raise ValueError(
"I expected all triangles to be cut along 1 or 2 edges. "
"I do not know how to proceed now.")
if return_intersection:
IS = Mesh(eltype='line2')
# Process the elements cutting one edge
#######################################
ncut1 = ncut==1
if ncut1.any():
prop1 = np.where(ncut1)[0]
fac1 = fac[ncut1]
ele1 = ele[ncut1]
cutedg1 = cutedg[fac1]
cut_edges = fac1[cutedg1]
# Identify the node numbers
# 0 : vertex on positive side
# 1 : vertex in plane
# 2 : new point dividing edge
# 3 : vertex on negative side
elems = np.column_stack([
ele1[p_pos[ele1]],
ele1[p_in[ele1]],
rev[cut_edges],
ele1[p_neg[ele1]],
])
if side in '+':
Sp += TriSurface(coords, elems[:, 0:3], prop=prop1)
if side in '-':
Sn += TriSurface(coords, elems[:, 1:4], prop=prop1)
# Process the elements cutting two edges
########################################
print("Cutting 2 edges")
ncut2 = ncut==2 # selector over whole range
print(ncut)
print(ncut2)
print(p_pos.sum(axis=-1)==2)
if ncut2.any():
prop2 = np.where(ncut2)[0]
fac2 = fac[ncut2]
ele2 = ele[ncut2]
pp2 = p_pos[ele2]
print("ele", ele2, pp2)
ncut2p = pp2.sum(axis=-1)==1 # selector over ncut2 elems
ncut2n = pp2.sum(axis=-1)==2
print(ncut2p, ncut2n)
if ncut2p.any():
prop1 = prop2[ncut2p]
fac1 = fac2[ncut2p]
ele1 = ele2[ncut2p]
print("ele1,1p", ele1)
cutedg1 = cutedg[fac1]
print(cutedg, fac1, cutedg1, fac1[cutedg1])
cut_edges = fac1[cutedg1].reshape(-1, 2)
corner = ele1[p_pos[ele1]]
quad = edg[cut_edges].reshape(-1, 4)
quad2 = quad[quad != corner.reshape(-1, 1)]
# Identify the node numbers
# 0 : vertex on positive side
# 1,2 : new points dividing edges
# 3,4 : vertices on negative side
elems = np.column_stack([
ele1[p_pos[ele1]],
rev[cut_edges],
quad2.reshape(-1, 2)
# ele1[p_neg[ele1]].reshape(-1,2),
])
if side in '+':
Sp += TriSurface(coords, elems[:, 0:3], prop=prop1)
if side in '-':
Sn += TriSurface(coords, elems[:, 1:4], prop=prop1)
Sn += TriSurface(coords, elems[:, 2:5], prop=prop1)
if ncut2n.any():
# print "# one vertex at negative side"
prop1 = np.where(ncut2n)[0]
fac1 = fac[ncut2n]
ele1 = ele[ncut2n]
cutedg1 = cutedg[fac1]
cut_edges = fac1[cutedg1].reshape(-1, 2)
# print cut_edges
corner = ele1[p_neg[ele1]]
# print corner
quad = edg[cut_edges].reshape(-1, 4)
# print quad
# print quad != corner.reshape(-1,1)
quad2 = quad[quad != corner.reshape(-1, 1)]
# print quad2
# 0 : vertex on negative side
# 1,2 : new points dividing edges
# 3,4 : vertices on positive side
elems = np.column_stack([
quad2.reshape(-1, 2),
# we can not use this, because order of the 2 vertices
# is importtant
# ele1[p_pos[ele1]].reshape(-1,2),
rev[cut_edges],
ele1[p_neg[ele1]],
])
# print elems
if side in '+':
Sp += TriSurface(coords, elems[:, 0:3], prop=prop1)
Sp += TriSurface(coords, elems[:, 1:4], prop=prop1)
if side in '-':
Sn += TriSurface(coords, elems[:, 2:5], prop=prop1)
return finalize(Sp, Sn, IS)
# Result
if side in '+':
Sp = Sp.compact() # .fixNormals()
if side in '-':
Sn = Sn.compact() # .fixNormals()
return _select_side(side, [Sp, Sn])
[docs] def cutWithPlane(self, *args, **kargs):
"""Cut a surface with a plane or a set of planes.
Cut the surface with one or more plane and returns either one side
or both. This uses a conversion to a 3-plex Formex to do the
cutting, and then converts the results back to TriSurfaces.
The parameters are the same as in :meth:`Formex.CutWithPlane`.
The returned surface(s) will have its normals fixed wherever possible.
"""
F = self.toFormex()
R = F.cutWithPlane(*args, **kargs)
if isinstance(R, list):
return [TriSurface(r).fixNormals() for r in R]
else:
return TriSurface(R).fixNormals()
def intersectionWithLines(self, p, p1, method='line', atol=1.e-5):
"""Compute the intersection points with a set of lines.
Parameters
----------
p: :term:`coords_like` (nlines,3)
A first point for each of the lines to intersect.
p1 :term:`coords_like` (nlines,3)
The second point defining the lines to intersect.
method: 'line' | 'segment' | ' ray'
Define whether the points ``p`` and ``p1`` define an infinitely
long line, a finite segment p-p1 or a half infinite ray (p->p1).
atol : float
Tolerance to be used in deciding whether an intersection point on
a border edge is inside the surface or not.
Returns
-------
X: Coords (nipts, 3)
A fused set of intersection points
ind: int array (nipts, 3)
An array identifying the related intersection points, lines and
triangle.
Notes
-----
A line laying in the plane of a triangle does not generate
intersections.
This method is faster than the similar function
:func:`geomtools.intersectionPointsLWT`.
Examples
-------
>>> T = Formex('3:.12.34').toSurface()
>>> L = Coords([[[0.,0.,0.], [0.,0.,1.]],
... [[0.5,0.5,0.5], [0.5,0.5,1.]],
... [[0.2,0.7,0.0], [0.2,0.7,1.]],
... ])
>>> P,X = T.intersectionWithLines(L[:,0,:], L[:,1,:])
>>> print(P)
[[ 0. 0. 0. ]
[ 0.5 0.5 0. ]
[ 0.2 0.7 0. ]]
>>> print(X)
[[0 0 0]
[0 0 1]
[1 1 0]
[1 1 1]
[2 2 1]]
"""
return geomtools.intersectLineWithTriangle(
self.coords[self.elems], p, p1, method=method, atol=atol)
[docs] def intersectionWithPlane(self, p, n, atol=0., sort='number'):
"""Return the intersection lines with plane (p,n).
Returns a plex-2 mesh with the line segments obtained by cutting
all triangles of the surface with the plane (p,n)
p is a point specified by 3 coordinates.
n is the normal vector to a plane, specified by 3 components.
The return value is a plex-2 Mesh where the line segments defining
the intersection are sorted to form continuous lines. The Mesh has
property numbers such that all segments forming a single continuous
part have the same property value.
By default the parts are assigned property numbers in decreasing
order of the number of line segments in the part. Setting the sort
argument to 'distance' will sort the parts according to increasing
distance from the point p.
The splitProp() method can be used to get a list of Meshes.
"""
n = np.asarray(n)
p = np.asarray(p)
# The vertices are classified based on their distance d from the plane:
# - inside: d = 0
# - up: d > 0
# - down: d < 0
# First, reduce the surface to the part intersecting with the plane:
# remove triangles with all up or all down vertices
d = self.distanceFromPlane(p, n)
d_ele = d[self.elems]
ele_all_up = (d_ele > 0.).all(axis=1)
ele_all_do = (d_ele < 0.).all(axis=1)
S = self.cselect(ele_all_up+ele_all_do, compact=False)
# If there is no intersection, we're done
if S.nelems() == 0:
return Mesh(Coords(), Connectivity(nplex=2, eltype='line2'))
Mparts = []
coords = S.coords
edg = S.getEdges()
fac = S.getElemEdges()
ele = S.elems
# No need to recompute distances, as clipping does not compact!
# Get the edges intersecting with the plane: 1 up and 1 down vertex
d_edg = d[edg]
edg_1_up = (d_edg > 0.).sum(axis=1) == 1
edg_1_do = (d_edg < 0.).sum(axis=1) == 1
w = edg_1_up * edg_1_do
ind = np.where(w)[0]
# compute the intersection points
if ind.size != 0:
rev = at.inverseUniqueIndex(ind)
M = Mesh(S.coords, edg[w])
x = geomtools.intersectionPointsSWP(
M.toFormex().coords, p, n, mode='pair',
return_all=True, atol=atol).reshape(-1, 3)
# For each triangle, compute the number of cutting edges
cut = w[fac]
ncut = cut.sum(axis=1)
# Split the triangles based on the number of inside vertices
d_ele = d[ele]
ins = d_ele == 0.
nins = ins.sum(axis=1)
ins0, ins1, ins2, ins3 = [np.where(nins==i)[0] for i in range(4)]
# No inside vertices -> 2 cutting edges
if ins0.size > 0:
cutedg = fac[ins0][cut[ins0]].reshape(-1, 2)
e0 = rev[cutedg]
Mparts.append(Mesh(x, e0, eltype='line2').compact())
# One inside vertex
if ins1.size > 0:
ncut1 = ncut[ins1]
cut10, cut11 = [np.where(ncut1==i)[0] for i in range(2)]
# 0 cutting edges: does not generate a line segment
# 1 cutting edge
if cut11.size != 0:
e11ins = ele[ins1][cut11][ins[ins1][cut11]].reshape(-1, 1)
cutedg = fac[ins1][cut11][cut[ins1][cut11]].reshape(-1, 1)
e11cut = rev[cutedg]
x11 = Coords.concatenate([coords, x], axis=0)
e11 = np.concatenate([e11ins, e11cut+len(coords)], axis=1)
Mparts.append(Mesh(x11, e11, eltype='line2').compact())
# Two inside vertices -> 0 cutting edges
if ins2.size > 0:
e2 = ele[ins2][ins[ins2]].reshape(-1, 2)
Mparts.append(Mesh(coords, e2, eltype='line2').compact())
# Three inside vertices -> 0 cutting edges
if ins3.size > 0:
insedg = fac[ins3].reshape(-1)
insedg.sort(axis=0)
double = insedg == np.roll(insedg, 1, axis=0)
insedg = np.setdiff1d(insedg, insedg[double])
if insedg.size != 0:
e3 = edg[insedg]
Mparts.append(Mesh(coords, e3, eltype='line2').compact())
# Done with getting the segments
if len(Mparts) == 0:
# No intersection: return empty mesh
return Mesh(Coords(), Connectivity(nplex=2, eltype='line2'))
else:
M = Mesh.concatenate(Mparts)
# Remove degenerate and duplicate elements
M = Mesh(M.coords, M.elems.removeDegenerate().removeDuplicate())
# Split in connected loops
parts = M.elems.chained()
prop = np.concatenate([[i]*part.nelems()
for i, part in enumerate(parts)])
elems = np.concatenate(parts, axis=0)
if sort == 'distance':
d = np.array([M.coords[part].distanceFromPoint(p).min()
for part in parts])
srt = np.argsort(d)
inv = at.inverseUniqueIndex(srt)
prop = inv[prop]
return Mesh(M.coords, elems, prop=prop)
[docs] def slice(self, dir=0, nplanes=20):
"""Intersect a surface with a sequence of planes.
A sequence of nplanes planes with normal dir is constructed
at equal distances spread over the bbox of the surface.
The return value is a list of intersectionWithPlane() return
values, i.e. a list of Meshes, one for every cutting plane.
In each Mesh the simply connected parts are identified by
property number.
"""
o = self.center()
if at.isInt(dir):
dir = at.unitVector(dir)
xmin, xmax = self.coords.directionalExtremes(dir, o)
P = Coords.interpolate(xmin, xmax, nplanes)
return [self.intersectionWithPlane(p, dir) for p in P]
################## Smooth a surface #############################
[docs] def smooth(self, method='lowpass', iterations=1, lambda_value=0.5,
neighbourhood=1, alpha=0.0, beta=0.2):
"""Smooth the surface.
Returns a TriSurface which is a smoothed version of the original.
Two smoothing methods are available: 'lowpass' and 'laplace'.
Parameters:
- `method`: 'lowpass' or 'laplace'
- `iterations`: int: number of iterations
- `lambda_value`: float: lambda value used in the filters
Extra parameters for 'lowpass' and 'laplace':
- `neighbourhood`: int: maximum number of edges followed in defining
the node neighbourhood
Extra parameters for 'laplace':
- `alpha`, `beta`: float: parameters for the laplace method.
Returns the smoothed TriSurface
"""
method = method.lower()
# find adjacency
adj = _adjacencyArrays(self.getEdges(), nsteps=neighbourhood)[1:]
adj = np.column_stack(adj)
# find interior vertices
bound_edges = self.borderEdgeNrs()
inter_vertex = np.resize(True, self.ncoords())
inter_vertex[np.unique(self.getEdges()[bound_edges])] = False
# calculate weights
w = np.ones(adj.shape, dtype=float)
w[adj<0] = 0.
val = (adj>=0).sum(-1).reshape(-1, 1)
w /= val
w = w.reshape(adj.shape[0], adj.shape[1], 1)
# recalculate vertices
if method == 'laplace':
xo = self.coords
x = self.coords.copy()
for step in range(iterations):
xn = x + lambda_value*(w*(x[adj]-x.reshape(-1, 1, 3))).sum(1)
xd = xn - (alpha*xo + (1-alpha)*x)
x[inter_vertex] = xn[inter_vertex] - (
beta*xd[inter_vertex]
+ (1-beta)*(w[inter_vertex]*xd[adj[inter_vertex]]).sum(1))
else: # default: lowpass
k = 0.1
mu_value = -lambda_value/(1-k*lambda_value)
x = self.coords.copy()
for step in range(iterations):
x[inter_vertex] = x[inter_vertex] + lambda_value*(
w[inter_vertex]*(x[adj[inter_vertex]]
- x[inter_vertex].reshape(-1, 1, 3))).sum(1)
x[inter_vertex] = x[inter_vertex] + mu_value*(
w[inter_vertex]*(x[adj[inter_vertex]]
- x[inter_vertex].reshape(-1, 1, 3))).sum(1)
return TriSurface(x, self.elems, prop=self.prop)
[docs] def smoothLowPass(self, iterations=2, lambda_value=0.5, neighbours=1):
"""Apply a low pass smoothing to the surface."""
return self.smooth('lowpass', iterations//2, lambda_value, neighbours)
[docs] def smoothLaplaceHC(self, iterations=2, lambda_value=0.5, alpha=0., beta=0.2):
"""Apply Laplace smoothing with shrinkage compensation to the surface."""
return self.smooth('laplace', iterations, lambda_value,
alpha=alpha, beta=beta)
[docs] def refine(self, max_edges=None, min_cost=None, method='gts'):
"""Refine the TriSurface.
Refining a TriSurface means increasing the number of triangles and
reducing their size, while keeping the changes to the modeled surface
minimal.
Construct a refined version of the surface.
This uses the external program `gtsrefine`. The surface
should be a closed orientable non-intersecting manifold.
Use the :meth:`check` method to find out.
Parameters:
- `max_edges`: int: stop the refining process if the number of
edges exceeds this value
- `min_cost`: float: stop the refining process if the cost of refining
an edge is smaller
- `log`: boolean: log the evolution of the cost
- `verbose`: boolean: print statistics about the surface
"""
if method == 'gts':
return self.gts_refine(max_edges, min_cost)
# THIS IS WORK IN PROGRESS
self.getElemEdges()
edglen = at.length(self.coords[self.edges[:, 1]]
- self.coords[self.edges[:, 0]])
print(edglen)
return self
[docs] def similarity(self, S):
"""Compute the similarity with another TriSurface.
Compute a quantitative measure of the similarity of the volumes
enclosed by two TriSurfaces. Both the calling and the passed
TriSurface should be closed manifolds (see :meth:`isClosedManifold`).
Returns a tuple (jaccard, dice, overlap).
If A and B are two closed manifolds, VA and VB are their respective
volumes, VC is the volume of the intersection of A and B, and VD is
the volume of the union of A and B, then the following similarity
measures are defined:
- jaccard coefficient: VC / VD
- dice: 2 * VC / (VA + VB)
- overlap: VC / min(VA,VB)
Both jaccard and dice range from 0 when the surfaces are completely
disjoint to 1 when the surfaces are identical. The overlap coefficient
becomes 1 when one of the surfaces is completely inside the other.
This method uses gts library to compute the intersection or union.
If that fails, nan values are returned.
"""
A, B = self, S
VA = A.volume()
VB = B.volume()
try:
VC = A.boolean(B, '*').volume()
VD = VA + VB - VC
except Exception:
try:
VD = A.boolean(B, '+').volume()
VC = VA + VB - VD
except Exception:
VC = VD = np.nan
dice = 2 * VC / (VA+VB)
overlap = VC / min([VA, VB])
jaccard = VC / VD
return jaccard, dice, overlap
###################################################################
## Methods using admesh/GTS
##############################
[docs] def fixNormals(self, outwards=True):
"""Fix the orientation of the normals.
Some surface operations may result in improperly oriented normals,
switching directions from one triangle to the adjacent one.
This method tries to reverse improperly oriented normals so that a
singly oriented surface is achieved.
.. note: In the current version, this uses the external program
`admesh`, so this should be installed on the machine.
If the surface is a (possibly non-orientable) manifold, the result
will be an orientable manifold.
If the surface is a closed manifold, the normals will be
oriented to the outside. This is done by computing the volume
inside the surface and reversing the normals if that turns out
to be negative.
Parameters:
- `outwards`: boolean: if True (default), a test is done whether
the surface is a closed manifold, and if so, the normals are
oriented outwards. Setting this value to False will skip this
test and the (possible) reversal of the normals.
"""
if self.nelems() == 0:
# Protect against impossible file handling for empty surfaces
return self
with utils.TempDir() as tmpdir:
tmp = tmpdir / 'surface.stl'
tmp1 = tmpdir / 'surface.off'
print(f"Writing temp file {tmp}")
self.write(tmp, 'stl')
print(f"Fixing surface while converting to OFF format {tmp1}")
fileread.stlConvert(tmp, tmp1)
print(f"Reading result from {tmp1}")
S = TriSurface.read(tmp1)
S.setProp(self.prop)
if outwards and S.isClosedManifold() and S.volume() < 0.:
S = S.reverse()
return S
[docs] def check(self, matched=True, verbose=False):
"""Check the surface using gtscheck.
Uses the external program `gtscheck` to check whether the surface
is an orientable, non self-intersecting manifold.
This is a necessary condition for using the `gts` methods:
split, coarsen, refine, boolean. Additionally, the surface should be
closed: this can be checked with :meth:`isClosedManifold`.
Parameters
----------
matched: bool
If True, self intersecting triangles are returned as element
indices of self. This is the default. If False, the self
intersecting triangles are returned as a separate TriSurface.
verbose: bool
If True, prints the statistics reported by the gtscheck
command.
Returns
-------
status: int
Return code from the checking program. One of the following:
- 0: the surface is an orientable, non self-intersecting manifold.
- 1: the created GTS file is invalid: this should normally not occur.
- 2: the surface is not an orientable manifold. This may be due to
misoriented normals. The :meth:`fixNormals` and :meth:`reverse`
methods may be used to help fixing the problem in such case.
- 3: the surface is an orientable manifold but is
self-intersecting. The self intersecting triangles are returned as
the second return value.
intersect: None | list of ints | TriSurface
None in case of a ``status`` 0, 1 or 2. For ``status`` value 3,
returns the self intersecting triangles as a list of element
numbers (if ``matched`` is True) or as a TriSurface (if ``matched``
is False).
"""
with utils.TempDir() as tmpdir:
tmp = tmpdir / 'surface.gts'
self.write(tmp, 'gts')
P = utils.system("gtscheck -v", stdin=tmp)
if verbose:
print(P.returncode)
if P.returncode == 0:
print("The surface is an orientable "
"non self-intersecting manifold")
res = None
elif P.returncode==2:
print("The surface is not an orientable manifold "
"(this may be due to badly oriented normals)")
res = None
elif P.returncode==3:
print("The surface is an orientable manifold "
"but is self-intersecting")
tmp = tmpdir / 'intersect.gts'
print(f"Writing temp file {tmp}")
with open(tmp, 'w') as fil:
fil.write(P.stdout)
res = TriSurface.read(tmp)
if matched:
res = self.matchCentroids(res)
else:
print("Status of gtscheck not understood")
res = None
return P.returncode, res
[docs] def split(self, base, verbose=False):
"""Split the surface using gtssplit.
Splits the surface into connected and manifold components.
This uses the external program `gtssplit`. The surface
should be a closed orientable non-intersecting manifold.
Use the :meth:`check` method to find out.
This method creates a series of files with given base name,
each file contains a single connected manifold.
"""
with utils.TempDir() as tmpdir:
tmp = tmpdir / 'surface.gts'
print(f"Writing surface to file {tmp}")
self.write(tmp, 'gts')
cmd = f"gtssplit -v {base}"
if verbose:
cmd += ' -v'
print(f"Splitting with command\n {cmd}")
P = utils.command(cmd, stdin=tmp, shell=True)
if P.returncode or verbose:
print(P.stdout)
#
# TODO: WE SHOULD READ THESE FILES BACK !!!
#
[docs] def coarsen(self, min_edges=None, max_cost=None,
mid_vertex=False, length_cost=False, max_fold=1.0,
volume_weight=0.5, boundary_weight=0.5, shape_weight=0.0,
progressive=False, log=False, verbose=False):
"""Coarsen a surface using gtscoarsen.
Construct a coarsened version of the surface.
This uses the external program `gtscoarsen`. The surface
should be a closed orientable non-intersecting manifold.
Use the :meth:`check` method to find out.
Parameters
----------
min_edges: int
Stop the coarsening process if the number of edges was to fall
below it.
max_cost: float
Stop the coarsening process if the cost of collapsing an edge is
larger than the specified value.
mid_vertex: bool
Use midvertex as replacement vertex instead of the default, which
is a volume optimized point.
length_cost: bool
Use length^2 as cost function instead of the default optimized
point cost.
max_fold: float
Maximum fold angle in degrees.
volume_weight: float
Weight used for volume optimization.
boundary_weight: float
Weight used for boundary optimization.
shape_weight: float
Weight used for shape optimization.
progressive: bool
If True, write progressive surface file.
log: bool
If Trye, log the evolution of the cost.
verbose: bool
If True, print statistics about the surface.
"""
if min_edges is None and max_cost is None:
min_edges = self.nedges() // 2
cmd = 'gtscoarsen'
if min_edges:
cmd += f" -n {min_edges}"
if max_cost:
cmd += f" -c {max_cost}"
if mid_vertex:
cmd += " -m"
if length_cost:
cmd += " -l"
if max_fold:
cmd += f" -f {max_fold}"
cmd += f" -w {volume_weight}"
cmd += f" -b {boundary_weight}"
cmd += f" -s {shape_weight}"
if progressive:
cmd += " -p"
if log:
cmd += " -L"
if verbose:
cmd += " -v"
with utils.TempDir() as tmpdir:
tmp = tmpdir / "surface.gts"
tmp1 = tmpdir / "surface1.gts"
print(f"Writing temp file {tmp}")
self.write(tmp, "gts")
print(f"Coarsening with command\n {cmd}")
P = utils.command(cmd, stdin=tmp, stdout=tmp1)
if P.returncode or verbose:
print(P.stdout)
print(f"Reading coarsened model from {tmp1}")
S = TriSurface.read(tmp1)
return S
[docs] def gts_refine(self, max_edges=None, min_cost=None, log=False, verbose=False):
"""Refine the TriSurface.
Refining a TriSurface means increasing the number of triangles and
reducing their size, while keeping the changes to the modeled surface
minimal.
This uses the external program `gtsrefine`. The surface
should be a closed orientable non-intersecting manifold.
Use the :meth:`check` method to find out.
Parameters
----------
max_edges: int
Stop the refining process if the number of edges exceeds this value.
min_cost: float
Stop the refining process if the cost of refining an edge is smaller.
log: bool
If True, log the evolution of the cost.
verbose: bool
If True, print statistics about the surface.
"""
if max_edges is None and min_cost is None:
max_edges = self.nedges() * 2
cmd = "gtsrefine"
if max_edges:
cmd += f" -n {max_edges}"
if min_cost:
cmd += f" -c {min_cost}"
if log:
cmd += " -L"
if verbose:
cmd += " -v"
with utils.TempDir() as tmpdir:
tmp = tmpdir / "surface.gts"
tmp1 = tmpdir / "surface1.gts"
print(f"Writing temp file {tmp}")
self.write(tmp, "gts")
print(f"Refining with command\n {cmd}")
P = utils.command(cmd, stdin=tmp, stdout=tmp1)
if P.returncode or verbose:
print(P.stdout)
print(f"Reading refined model from {tmp1}")
S = TriSurface.read(tmp1)
return S
[docs] def gts_smooth(self, iterations=1, lambda_value=0.5, verbose=False):
"""Smooth the surface using gtssmooth.
Smooth a surface by applying iterations of a Laplacian filter.
This uses the external program `gtssmooth`. The surface
should be a closed orientable non-intersecting manifold.
Use the :meth:`check` method to find out.
Parameters
----------
lambda_value: float
Laplacian filter parameter.
iterations: int
Number of iterations.
verbose: bool
If True, print statistics about the surface.
See Also
--------
smoothLowPass
smoothLaplaceHC
"""
cmd = "gtssmooth"
# if fold_smoothing:
# cmd += f" -f {fold_smoothing}"
cmd += f" {lambda_value} {iterations}"
if verbose:
cmd += " -v"
with utils.TempDir() as tmpdir:
tmp = tmpdir / "surface.gts"
tmp1 = tmpdir / "surface1.gts"
print(f"Writing temp file {tmp}")
self.write(tmp, "gts")
print(f"Smoothing with command\n {cmd}")
P = utils.command(cmd, stdin=tmp, stdout=tmp1)
if P.returncode or verbose:
print(P.stdout)
print(f"Reading smoothed model from {tmp1}")
S = TriSurface.read(tmp1)
return S
[docs] def gts_set(self, surf, op, prop=[1, 1, 2, 2], check=False, verbose=False):
"""Perform a boolean operation with another surface.
Boolean operations between surfaces are a basic operation in
free surface modeling. Both surfaces should be closed orientable
non-intersecting manifolds. Use the :meth:`check` method to find out.
Following is a list of defined operations, where surface 1 relates to
`self` and surface 2 to the `surf` argument. For simplicity, the
operations are identified by a short string.
Basic operations:
'i' : the part of surface 1 inside surface 2
'o' : the part of surface 1 outside surface 2
'2i' : the part of surface 2 inside surface 1
'2o' : the part od surface 2 outside surface 1
These surface are manifolds, but may be closed or not.
From these basic parts, the following set operation can be constructed.
These are mathematical set operation on the volumes in side the
surfaces, and result in closed surfaces:
'+' : the union of surfaces 1 and 2
'*' : the intersection of surfaces 1 and 2
'-' : the difference of surface 1 minus surface 2
'2-' : the difference of surface 2 minus surface 1
'^' : the symmetric difference of the surfaces (equal to '-' + '2-')
Parameters
----------
surf: TriSurface
Another TriSurface that is a closed manifold surface.
op: str or list of str
The operation(s) to perform: one of the operations specified
above, or a list of such operations. A special value 'a' will
return the full list of 9 surfaces in the above order.
prop: list of int
A list of 4 integer values that will be set as props on the
four base surfaces, to facilitate identification of the parts
of the result(s). The default value will give prop values 1 or
2 depending on the original surface the parts belonged to.
Specifying None or an empty list will return surfaces without
props.
check: bool
If True, a check is done that the surfaces are not self-intersecting;
if one of them is, the set of self-intersecting faces is written
(as a GtsSurface) on standard output
verbose: bool
If True, print statistics about the surface.
Returns
-------
:class:`TriSurface` or list thereof
A single manifold surface, or a list of such surfaces, corresponding
to the specified oppetaion(s). The base operation may not be closed.
The set operations always are closed.
Note
----
This method uses the external command 'gtsset' and will not run if
it is not installed (available from pyformex/extras).
"""
from pyformex.plugins.pyformex_gts import gtsset
base = gtsset(self, surf, op='a', filt='', ext='.gts',
check=check, verbose=verbose)
if not base:
raise ValueError("Error computing the base surfaces.")
# Function to compute the requested surface(s)
if prop:
base['s1in2'].setProp(prop[0])
base['s1out2'].setProp(prop[1])
base['s2in1'].setProp(prop[2])
base['s2out1'].setProp(prop[3])
def getres(o):
if o == 'i':
return base['s1in2']
elif o == 'o':
return base['s1out2']
elif o == '2i':
return base['s2in1']
elif o == '2o':
return base['s2out1']
elif o == '+':
return base['s1out2'] + base['s2out1']
elif o == '*':
return base['s1in2'] + base['s2in1']
elif o == '-':
return base['s1out2'] + base['s2in1'].reverse()
elif o == '2-':
return base['s2out1'] + base['s1in2'].reverse()
elif o == '^':
return (base['s1out2'] + base['s2out1']
+ (base['s1in2'] + base['s2in1']).reverse())
if op == 'a':
op = ['i', 'o', '2i', '2o', '+', '*', '-', '2-', '^']
if utils.isString(op):
return getres(op)
else:
return [getres(o) for o in op]
[docs] def boolean(self, surf, op, check=False, verbose=False):
"""Perform a boolean operation with another surface.
Boolean operations between surfaces are a basic operation in
free surface modeling. Both surfaces should be closed orientable
non-intersecting manifolds.
Use the :meth:`check` method to find out.
The boolean operations are set operations on the enclosed volumes:
union('+'), difference('-') or intersection('*').
Parameters
----------
surf: TriSurface
Another TriSurface that is a closed manifold surface.
op: '+', '-' or '*'
The boolean operation to perform: union('+'), difference('-')
or intersection('*').
check: bool
If True, a check is done that the surfaces are not self-intersecting;
if one of them is, the set of self-intersecting faces is written
(as a GtsSurface) on standard output
verbose: bool
If True, print statistics about the surface.
Returns
-------
TriSurface
A closed manifold TriSurface that is the volume union, difference or
intersection of self with surf.
Note
----
This method uses the external command 'gtsset' and will not run if
it is not installed (available from pyformex/extras).
"""
from pyformex.plugins.pyformex_gts import gtsset
return gtsset(self, surf, op, filt='', ext='.gts',
check=check, verbose=verbose).fuse().compact()
[docs] def intersection(self, surf, check=False, verbose=False):
"""Return the intersection curve(s) of two surfaces.
Boolean operations between surfaces are a basic operation in
free surface modeling. Both surfaces should be closed orientable
non-intersecting manifolds.
Use the :meth:`check` method to find out.
Parameters
----------
surf: TriSurface
A closed manifold surface.
check: bool, optional
If True, a check is made that the surfaces are not self-intersecting;
if one of them is, the set of self-intersecting faces is written
(as a GtsSurface) on standard output
verbose: bool, optional
If True, statistics about the surface are printed on stdout.
Returns
-------
Mesh
A Mesh with eltype Line2 holding all the line segments of the
intersection curve(s).
"""
from pyformex.plugins.pyformex_gts import gtsset
return gtsset(self, surf, op='*', ext='.list', curve=True,
check=check, verbose=verbose)
[docs] def inside(self, pts, method='gts', tol='auto', multi=False):
"""Test which of the points pts are inside the surface.
Parameters
----------
pts: :term_`coords_like`
The points to check agains the surface.
method`: str
Method to be used for the detection. Depending on
the software you have installed the following are possible:
- 'gts': provided by pyformex-extra (default)
- 'vtk': provided by python-vtk (slower)
tol: float
Tolerance on equality of floating point values.
Returns
-------
int array
The indices of the points that are inside the surface.
The indices refer to the onedimensional list
of points as obtained from Coords(pts).points().
"""
pts = Coords(pts).points()
if method == 'gts':
from pyformex.plugins.pyformex_gts import inside
return inside(self, pts, tol, multi=multi)
elif method == 'vtk':
from pyformex.plugins.vtk_itf import inside
return inside(self, pts, tol)
[docs] def outside(self, pts, **kargs):
"""Returns the points outside the surface.
This is the complement of :meth:`inside`. See there for
parameters and return value.
"""
return at.complement(self.inside(pts, **kargs), len(pts))
[docs] def voxelize(self, n, bbox=0.01, return_formex=False):
"""Voxelize the volume inside a closed surface.
Parameters
----------
n: int or (int, int, int)
Resolution, i.e. number of voxel cells to use along the three axes.
If a single int is specified, the number of cells will be adapted
according to the surface's :meth:`sizes` (as the voxel cells are
always cubes). The specified number of voxels will be used along the
largest direction.
bbox: float or (point,point)
Defines the bounding box of the volume that needs to be voxelized.
A float specifies a relative amount to add to the surface's bounding
box. Note that this defines the bounding box of the centers of the
voxels.
return_formex: bool
If True, also returns a Formex with the centers of the voxels.
Returns
-------
voxels: int array (nz,ny,nx)
The array has a value 1 for the voxels whose center is inside the
surface, else 0.
centers: Formex
A plex-1 Formex with the centers of the voxels, and property values
0 or 1 if the point is respectively outside or inside the surface.
The voxel cell ordering in the Formex is z-direction first, then y,
then x.
Notes
-----
See also example Voxelize, for saving the voxel values in a stack
of binary images.
"""
if not self.isClosedManifold():
raise ValueError("The surface is non a closed manifold")
from pyformex import simple
if at.isFloat(bbox):
a, b = 1.0+bbox, bbox
bbox = self.bbox()
bbox = [a*bbox[0]-b*bbox[1], a*bbox[1]-b*bbox[0]]
bbox = at.checkArray(bbox, shape=(2, 3), kind='f')
if at.isInt(n):
sz = bbox[1]-bbox[0]
step = sz.max() / (n-1)
n = np.ceil(sz / step).astype(at.Int)
n = at.checkArray(n, shape=(3,), kind='i')
X = simple.regularGrid(bbox[0], bbox[0]+n*step, n, swapaxes=True)
ind = self.inside(X)
vox = np.zeros(n+1, dtype=np.uint8)
vox.ravel()[ind] = 1
if return_formex:
P = Formex(X.reshape(-1, 3))
P.setProp(vox.ravel())
return vox, P
return vox
[docs] def tetgen(self, quality=2.0, volume=None, filename=None):
"""Create a tetrahedral mesh inside the surface.
This uses :func:`~plugins.tetgen.tetMesh` to generate a quality
tetrahedral mesh inside the surface. The surface should be a closed
manifold.
Parameters
----------
quality: float
The quality of the output tetrahedral mesh. The value is a
constraint on the circumradius-to-shortest-edge ratio. The
default (2.0) already provides a high quality mesh. Providing
a larger value will reduce quality but increase speed. With
quality=None, no quality constraint will be imposed.
volume: float, optional
If provided, applies a maximum tetrahedron volume constraint.
filename: :term:`path_like`
Specifies where the intermediate files will be stored.
The default will use a temporary directory which will be destroyed
after return. If the path of an existing directory is provided,
the files will be stored in that directory with a name 'surface.off'
for the original surface model and files 'surface.1.*' for the
generated tetrahedral model (in tetgen format).
If the path does not exist or is an existing file, the parent
directory should exist and files are stored with the given
file name as base. Existing files will be silently overwritten.
Returns
-------
Mesh
A tetrahedral Mesh (eltype='tet4') filling the input surface,
provided the :func:`~plugins.tetgen.tetMesh` function finished
successfully.
"""
from pyformex.plugins import tetgen
if filename is None:
tmpdir = utils.TempDir()
fn = Path(tmpdir.name) / 'surface.off'
else:
filename = Path(filename)
if filename.exists() and filename.is_dir():
fn = filename / 'surface.off'
elif filename.parent.exists():
fn = filename.with_suffix('.off')
else:
raise ValueError(f"Invalid filename specified: {filename}")
self.write(fn)
res = tetgen.tetMesh(fn, quality, volume)
return res['tetgen.ele']
# Set TriSurface methods defined elsewhere
from pyformex.plugins.webgl import surface2webgl
webgl = surface2webgl
@utils.deprecated_by('TriSurface.facetArea', 'TriSurface.areas')
def facetArea(self):
return self.areas()
############################################################################
[docs]def fillBorder(border, method='radial', dir=None):
"""Create a triangulated surface inside a given closed polygonal line.
Parameters
----------
border: :class:`PolyLine`, :class:`Mesh` or :class:`Coords`
A closed polygonal line that forms the border of the triangulated
surface to be created. The polygon does not have to be planar.
The line can be provided as one of the following:
- a closed PolyLine,
- a 2-plex Mesh, with a Connectivity table such that the elements
in the specified order form a closed polyline,
- a simple Coords holding the subsequent vertices of the polygonal
border line.
method: str
Specifies the algorithm to be used to fill the polygon. Currently
available are:
- 'radial': this method adds a central point and connects all border
segments with the center to create triangles.
- 'border': this method creates subsequent triangles by connecting the
endpoints of two consecutive border segments and thus works its way
inwards until the hole is closed. Triangles are created at the line
segments that form the smallest angle.
See also Notes below.
Returns
-------
TriSurface
A TriSurface filling the hole inside the border.
Notes
-----
The 'radial' method produces nice results if the border is relative smooth,
nearly convex and nearly planar. It adds an extra point though, which may
be unwanted. On irregular 3D borders there is a high change that the
resulting TriSurface contains intersecting triangles.
The 'border' method is slower on large borders, does not introduce any
new point and has a better chance of avoiding intersecting triangles
on irregular 3D borders.
The resulting surface can be checked for intersecting triangles with the
:meth:`check` method.
Because the 'border' does not create any new points, the returned surface
can use the same point coordinate array as the input object.
"""
from pyformex.plugins.curve import PolyLine
if isinstance(border, Mesh) and border.nplex()==2:
if method == 'radial':
border = border.compact()
coords = border.coords
elems = border.elems[:, 0]
elif isinstance(border, PolyLine):
coords = border.coords
elems = None
elif isinstance(border, Coords):
coords = border.reshape(-1, 3)
elems = None
else:
raise ValueError("Expected a 2-plex Mesh, "
"a PolyLine or a Coords array as first argument")
if elems is None:
elems = np.arange(coords.shape[0])
n = elems.shape[0]
if n < 3:
raise ValueError("Expected at least 3 points.")
if method == 'radial':
coords = Coords.concatenate([coords, coords.center()])
elems = np.column_stack([elems, np.roll(elems, -1),
n*np.ones(elems.shape[0], dtype=at.Int)])
elif method == 'border':
# creating elems array at once (more efficient than appending)
tri = -np.ones((n-2, 3), dtype=at.Int)
# compute all internal angles
x = coords[elems]
e = np.arange(n)
v = np.roll(x, -1, axis=0) - x
v = at.normalize(v)
c = at.vectorPairCosAngle(np.roll(v, 1, axis=0), v)
# loop in order of smallest angles
itri = 0
while n > 3:
# find minimal angle
j = c.argmin()
i = (j - 1) % n
k = (j + 1) % n
tri[itri] = [e[i], e[j], e[k]]
# remove the point j of triangle i,j,k
# recompute adjacent angles of edge i,k
ii = (i-1) % n
v1 = at.normalize([v[e[ii]], x[e[k]] - x[e[i]]])
v2 = at.normalize([x[e[k]] - x[e[i]], v[e[k]]])
cnew = at.vectorPairCosAngle(v1, v2)
c = np.roll(np.concatenate([cnew, np.roll(c, 1-j)[3:]]), j-1)
e = np.roll(np.roll(e, -j)[1:], j)
n -= 1
itri += 1
tri[itri] = e
elems = elems[tri]
elif method == 'planar':
from pyformex.plugins import polygon
x = coords[elems]
e = np.arange(x.shape[0])
if dir is None:
dir = geomtools.smallestDirection(x)
X, C, A, a = polygon.projected(x, dir)
P = polygon.Polygon(X)
if P.area() < 0.0:
P = P.reverse()
e = at.reverseAxis(e, 0)
S = P.fill()
e = e[S.elems]
elems = elems[e]
else:
raise ValueError("Strategy should be either 'radial', 'border' or 'planar'")
return TriSurface(coords, elems)
##########################################################################
####### Deprecated functions #######################
@utils.deprecated_by('trisurface.Sphere', 'simple.sphere')
def Sphere(level=4, verbose=False, filename=None):
from pyformex import simple
return simple.sphere(ndiv=2**(level-1))
@utils.deprecated_by('trisurface.intersectSurfaceWithLines',
'trisurface.intersectionWithLines')
def intersectSurfaceWithLines(ts, qli, mli, atol=1.e-5):
pass
# this function has been removed. You can get the same results using:
# P,X = intersectionWithLines(self, q, m, q2,atol=1.e-5)
# P,L,T = P[X[:, 0]], X[:, 1], X[:, 2]
@utils.deprecated_by('trisurface.intersectSurfaceWithSegments',
'trisurface.intersectionWithLines')
def intersectSurfaceWithSegments(s1, segm, atol=1.e-5):
pass
# this function has been removed. You can get the same results using:
# P,X = intersectionWithLines(self, q=segm[:, 0], m=None, q2=segm[:, 1],
# atol=1.e-5)
# P,L,T = P[X[:, 0]], X[:, 1], X[:, 2]
def read_vtk_surface(fn):
try:
from pyformex.plugins import vtk_itf
except ImportError:
utils.warn("I could not import VTK. This probably means that Python"
"VTK bindings are not installed on your machine.")
return vtk_itf.read_vtk_surface(fn)
# End
