#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""This module holds a tools to clip meshes against a plane"""
from .mesh import *
[docs]class MeshClipper(object):
"""A class to perform mesh clipping operations.
Parameters
----------
source_mesh : Mesh
The mesh to be clipped.
plane : Plane, optional
The clipping plane. By default, the plane is the Oxy plane.
vicinity_tol : float, optional
The absolute tolerance to consider en vertex is on the plane. Default is 1e-3.
assert_closed_boundaries : bool, optional
False by default. When True, the mesh clipper will raise an exception if intersections with the clipping
plane are not closed. It may be caused by a non-watertight mesh.
verbose : bool, optional
False by default. If True, some messages on operations that are handled are printed.
"""
def __init__(self, source_mesh, plane=Plane(), vicinity_tol=1e-3, assert_closed_boundaries=False, verbose=False):
self._source_mesh = source_mesh
self._plane = plane
self._vicinity_tol = vicinity_tol
self._assert_closed_boundaries = assert_closed_boundaries
self._verbose = verbose
self.__internals__ = dict()
self._update()
@property
def verbose(self):
"""Get the current verbosity"""
return self._verbose
[docs] def verbose_on(self):
"""Switches ON the verbosity of the clipper."""
self._verbose = True
[docs] def verbose_off(self):
"""Switches OFF the verbosity of the clipper."""
self._verbose = False
@property
def assert_closed_boundaries(self):
"""Do we assert the boundaries have to be closed"""
return self._assert_closed_boundaries
[docs] def assert_closed_boundaries_on(self):
"""Switches ON the flag for closed boundary assertion while clipping."""
self._assert_closed_boundaries = True
[docs] def assert_closed_boundaries_off(self):
"""Switches OFF the flag for closed boundary assertion while clipping."""
self._assert_closed_boundaries = False
return
@property
def vicinity_tol(self):
"""Vicinity tolerance.
It tells if a point is close enough to the plane to consider it lies on the plane
"""
return self._vicinity_tol
@vicinity_tol.setter
def vicinity_tol(self, value):
"""Set the vicinity tolerance that tells that a vertex is located on the plane."""
self.__internals__.clear()
self._vicinity_tol = float(value)
self._update()
@property
def source_mesh(self):
"""The mesh we work with"""
return self._source_mesh
@property
def plane(self):
"""The clipping plane"""
return self._plane
@plane.setter
def plane(self, value):
"""Changes the clipping plane."""
self.__internals__.clear()
self._plane = value
self._update()
def _update(self):
"""Updates the clipper"""
self._vertices_positions_wrt_plane()
self._partition_mesh()
self._clip()
def _vertices_positions_wrt_plane(self):
"""
Classifies vertices with respect to the clipping plane
"""
vertices_distances = self._plane.get_point_dist_wrt_plane(self._source_mesh.vertices)
vertices_positions = {'vertices_distances': vertices_distances,
'vertices_above_mask': vertices_distances > self._vicinity_tol,
'vertices_on_mask': np.fabs(vertices_distances) < self._vicinity_tol,
'vertices_below_mask': vertices_distances < -self._vicinity_tol
}
if np.all(vertices_positions['vertices_above_mask']):
# Cannot clip as this mesh is totally above or below the clipping plane. Throwing an error
raise RuntimeError("Mesh does not crosses the clipping plane", "above")
if np.all(vertices_positions['vertices_below_mask']):
raise RuntimeError("Mesh does not crosses the clipping plane", "below")
self.__internals__.update(vertices_positions)
def _partition_mesh(self):
"""Partitions the mesh in 3 with respect to the plane
* upper_mesh: part entirely above the clipping plane
* crown_mesh: part intersecting the clipping plane
* lower_mesh: part entirely under the clipping plane
"""
vertices_distances = self.__internals__['vertices_distances']
vertices_above_mask = self.__internals__['vertices_above_mask']
vertices_on_mask = self.__internals__['vertices_on_mask']
vertices_below_mask = self.__internals__['vertices_below_mask']
source_mesh_faces = self.source_mesh.faces
nb_vertices_above = vertices_above_mask[source_mesh_faces].sum(axis=1)
nb_vertices_below = vertices_below_mask[source_mesh_faces].sum(axis=1)
# Simple criteria ensuring that _faces are totally above or below the plane (4 _vertices at the same side)
# Works for both triangles and quadrangles
above_faces_mask = nb_vertices_above == 4
below_faces_mask = nb_vertices_below == 4
crown_faces_mask = np.logical_and(np.logical_not(above_faces_mask), np.logical_not(below_faces_mask))
above_faces_ids = np.where(above_faces_mask)[0]
below_faces_ids = np.where(below_faces_mask)[0]
crown_faces_ids = np.where(crown_faces_mask)[0]
partition = dict()
def generate_mesh(faces_ids, key):
new_mesh, ids = self._source_mesh.extract_faces(faces_ids, return_index=True)
new_mesh.name = key
partition['_'.join((key, 'vertices_ids'))] = ids
partition['_'.join((key, 'vertices_distances'))] = vertices_distances[ids]
partition['_'.join((key, 'above_vertices_mask'))] = vertices_above_mask[ids]
partition['_'.join((key, 'on_vertices_mask'))] = vertices_on_mask[ids]
partition['_'.join((key, 'below_vertices_mask'))] = vertices_below_mask[ids]
partition[key] = new_mesh
# Generating partition meshes
generate_mesh(above_faces_ids, 'upper_mesh')
generate_mesh(crown_faces_ids, 'crown_mesh')
generate_mesh(below_faces_ids, 'lower_mesh')
partition['above_faces_ids'] = above_faces_ids
partition['crown_faces_ids'] = crown_faces_ids
partition['below_faces_ids'] = below_faces_ids
self.__internals__.update(partition)
@property
def lower_mesh(self):
"""A new mesh composed of the faces that entirely lie under the clipping plane.
Returns
-------
Mesh
"""
return self.__internals__['lower_mesh']
@property
def crown_mesh(self):
"""A new mesh only having the faces that cut the plane
Returns
-------
Mesh
"""
return self.__internals__['crown_mesh']
@property
def upper_mesh(self):
"""A new mesh only having the faces lying entirely up the plane
Returns
-------
Mesh
"""
return self.__internals__['upper_mesh']
@property
def clipped_crown_mesh(self):
"""A new mesh that is obtained by clipping the crown_mesh.
Returns
-------
Mesh
"""
return self.__internals__['clipped_crown_mesh']
@property
def clipped_mesh(self):
"""The resulting clipped mesh"""
return self.__internals__['clipped_mesh']
@property
def closed_polygons(self):
"""Returns the list of closed boundary polygons obtained after clipping.
This is a list of lists. The enclosed lists are ordered IDs list that form a closed polygon, described in the
counter-clockwise order with respect to the mesh (oriented following the clipping plane's normal).
Returns
-------
list
Warnings
--------
* The first vertex is repeated at the end of the list. By definition, these polygons are lying on the clipping
plane.
* Vertices IDs are corresponding to the IDs of the clipped_crown_mesh, not those of the clipped_mesh.
"""
return self.__internals__['closed_polygons']
@property
def closed_polygons_vertices(self):
"""Returns the list of closed boundary polygons obtained after clipping.
This is a list of lists. The enclosed lists are ordered vertices coordinates of the closed polygons. By
definition, these polygons are lying on the clipping plane.
Returns
-------
list
"""
polygons = self.__internals__['closed_polygons']
closed_polygons_vertices = []
# TODO: voir si on ne peut pas directement indicer par polygons sans boucle for
for polygon in polygons:
closed_polygons_vertices.append(self.clipped_crown_mesh.vertices[polygon])
return closed_polygons_vertices
@property
def nb_closed_polygons(self):
"""The number of closed polygons obtained after clipping
Returns
-------
int
"""
return len(self.__internals__['closed_polygons'])
@property
def open_lines(self):
"""Returns a list of open boundary lines obtained after clipping.
This is a list of lists. The enclosed lists are ordered vertices IDs.
Returns
-------
list
Warning
-------
* The vertices IDs correspond to the IDs of clipped_crown_mesh, not clipped_mesh.
"""
return self.__internals__['open_lines']
@property
def open_lines_vertices(self):
"""Returns a list of open boundary lines obtained after clipping.
This is a list of lists. The enclosed lists are ordered vertices IDs.
Returns
-------
list
"""
lines = self.__internals__['open_lines']
lines_vertices = []
# TODO: voir si on ne peut pas directement indicer par polygons sans boucle for
for line in lines:
lines_vertices.append(self.clipped_crown_mesh.vertices[line])
return lines_vertices
@property
def nb_open_lines(self):
"""The number of open lines obtained after clipping.
Returns
-------
int
"""
return len(self.__internals__['open_lines'])
def _clip_crown_by_plane(self):
"""Performs the clipping operation of the crown_mesh and determines the obtained boundaries.
This is the heart method of the class.
"""
crown_mesh = self.crown_mesh
vertices = crown_mesh.vertices
vertices_on_mask = self.__internals__['crown_mesh_on_vertices_mask']
# TODO: Vertices pre-projection to be done here !!!
# vertices_on = partition['vertices_on']
# _vertices[vertices_on] = plane.orthogonal_projection_on_plane(_vertices[vertices_on])
# pos[vertices_on] = 0.
vertices_above_mask = self.__internals__['crown_mesh_above_vertices_mask']
vertices_below_mask = self.__internals__['crown_mesh_below_vertices_mask']
vertices_distances = self.__internals__['crown_mesh_vertices_distances']
# Init
crown_faces = list()
direct_boundary_edges = dict()
inv_boundary_edges = dict()
intersections = list()
index = crown_mesh.nb_vertices
for face_id in range(crown_mesh.nb_faces):
face = crown_mesh.get_face(face_id)
# # Determining the type of face clipping
v_above_face = np.where(vertices_above_mask[face])[0]
v_on_face = np.where(vertices_on_mask[face])[0]
v_below_face = np.where(vertices_below_mask[face])[0]
nb_above = len(v_above_face)
nb_on = len(v_on_face)
nb_below = len(v_below_face)
face_type = str(nb_above) + str(nb_on) + str(nb_below)
if face_type == '202': # Done
# 0*-----*3
# | |
# ----o-----o----
# | |
# 1*-----*2
if v_above_face[1] == v_above_face[0] + 1:
face = np.roll(face, -v_above_face[1])
p0, p1, p2, p3 = vertices[face]
ileft = self._plane.get_edge_intersection(p0, p1)
iright = self._plane.get_edge_intersection(p2, p3)
intersections += [ileft, iright]
boundary_edge = [index, index + 1]
crown_faces.append([index, face[1], face[2], index + 1])
index += 2
elif face_type == '301': # Done
# *2
# / \
# / \
# / \
# 3* *1
# \ /
# ---o---o---
# \ /
# *0
face = np.roll(face, -v_below_face[0])
p0, p1, p3 = vertices[face[[0, 1, 3]]]
ileft = self._plane.get_edge_intersection(p0, p3)
iright = self._plane.get_edge_intersection(p0, p1)
intersections += [ileft, iright]
boundary_edge = [index, index + 1]
crown_faces.append([index, face[0], index + 1, index])
index += 2
elif face_type == '103': # Done
# *0
# / \
# ---o---o---
# / \
# 1* - - - *3
# \ /
# \ /
# \ /
# *2
face = np.roll(face, -v_above_face[0])
p0, p1, p3 = vertices[face[[0, 1, 3]]]
ileft = self._plane.get_edge_intersection(p0, p1)
iright = self._plane.get_edge_intersection(p0, p3)
intersections += [ileft, iright]
boundary_edge = [index, index + 1]
crown_faces.append([index, face[1], face[3], index + 1])
crown_faces.append([face[1], face[2], face[3], face[1]])
index += 2
elif face_type == '102': # Done
# *O
# / \
# ---o---o---
# / \
# 1*-------*2
face = np.roll(face, -v_above_face[0])
p0, p1, p2 = vertices[face]
ileft = self._plane.get_edge_intersection(p0, p1)
iright = self._plane.get_edge_intersection(p0, p2)
intersections += [ileft, iright]
boundary_edge = [index, index + 1]
crown_faces.append([index, face[1], face[2], index + 1])
index += 2
elif face_type == '201': # done
# 2*-------*1
# \ /
# ---o---o---
# \ /
# *0
face = np.roll(face, -v_below_face[0])
p0, p1, p2 = vertices[face]
ileft = self._plane.get_edge_intersection(p0, p2)
iright = self._plane.get_edge_intersection(p0, p1)
intersections += [ileft, iright]
boundary_edge = [index, index + 1]
crown_faces.append([index, face[0], index + 1, index])
index += 2
elif face_type == '211': # Done
# *3 *1
# / \ / \
# / *2 or 2* \
# 0/ / \ \0
# ---*---o--- ---o---*---
# \ / \ /
# *1 *3
#
face = np.roll(face, -v_on_face[0])
if vertices_distances[face[1]] < 0.:
p1, p2 = vertices[face[[1, 2]]]
iright = self._plane.get_edge_intersection(p1, p2)
intersections.append(iright)
boundary_edge = [face[0], index]
crown_faces.append([face[0], face[1], index, face[0]])
else:
p2, p3 = vertices[face[[2, 3]]]
ileft = self._plane.get_edge_intersection(p2, p3)
intersections.append(ileft)
boundary_edge = [index, face[0]]
crown_faces.append([index, face[3], face[0], index])
index += 1
elif face_type == '112': # Done
# *3 *1
# / \ / \
# ---*---o--- or ---o---*---
# 0\ \ / /0
# \ *2 2* /
# \ / \ /
# *1 *3
face = np.roll(face, -v_on_face[0])
if vertices_distances[face[1]] < 0.:
p2, p3 = vertices[face[[2, 3]]]
iright = self._plane.get_edge_intersection(p2, p3)
intersections.append(iright)
boundary_edge = [face[0], index]
crown_faces.append([face[0], face[1], face[2], index])
else:
p1, p2 = vertices[face[[1, 2]]]
ileft = self._plane.get_edge_intersection(p1, p2)
intersections.append(ileft)
boundary_edge = [index, face[0]]
crown_faces.append([index, face[2], face[3], face[0]])
index += 1
elif face_type == '013': # Done
# -----*-----
# / \
# / \
# * *
# \ /
# \ /
# *
boundary_edge = None
crown_faces.append(list(face))
elif face_type == '210' or face_type == '310': # Done
# *-------* *
# \ 210 / / \ 310
# \ / * *
# \ / \ /
# ----*---- ----*----
boundary_edge = None
elif face_type == '111': # Done
# *2 *1
# /| |\
# / | | \
# ---*--o--- or ---o--*---
# 0\ | | /0
# \| |/
# *1 *2
face = np.roll(face, -v_on_face[0])
p1, p2 = vertices[face[[1, 2]]]
if vertices_distances[face[1]] < 0.:
iright = self._plane.get_edge_intersection(p1, p2)
intersections.append(iright)
boundary_edge = [face[0], index]
crown_faces.append([face[0], face[1], index, face[0]])
else:
ileft = self._plane.get_edge_intersection(p1, p2)
intersections.append(ileft)
boundary_edge = [index, face[0]]
crown_faces.append([index, face[2], face[0], index])
index += 1
elif face_type == '120': # Done
# *O
# / \
# / \
# 1/ \2
# ----*-------*----
# face = np.roll(face, -v_above_face[0])
# boundary_edge = [face[1], face[2]]
# FIXME: quick fix here : robust ?
boundary_edge = None
elif face_type == '021': # Done
# ----*-------*----
# 2\ /1
# \ /
# \ /
# *0
face = np.roll(face, -v_below_face[0])
boundary_edge = [face[2], face[1]]
face = list(face)
face.append(face[0])
crown_faces.append(face)
elif face_type == '022':
# ----*-----*----
# 0| |3
# | |
# 1*-----*2
if v_on_face[1] == v_on_face[0] + 1:
face = np.roll(face, -v_on_face[1])
boundary_edge = [face[0], face[3]]
crown_faces.append(list(face))
elif face_type == '012': # Done
# ------*------
# / \
# / \
# / \
# *-------*
boundary_edge = None
face = list(face)
face.append(face[0])
crown_faces.append(face)
elif face_type == '220': # Done
# 0*-----*3
# | |
# 1| |2
# ----*-----*----
# if v_above_face[1] == v_above_face[0] + 1:
# face = np.roll(face, -v_above_face[1])
# boundary_edge = [face[1], face[2]]
# FIXME: quick fix here : robust ?
boundary_edge = None
elif face_type == '121': # Done
# *0
# / \
# / \
# ---*-----*---
# 1\ /3
# \ /
# *2
face = np.roll(face, -v_above_face[0])
boundary_edge = [face[1], face[3]]
crown_faces.append([face[1], face[2], face[3], face[1]])
elif face_type == '300' or face_type == '400':
# * *-----*
# / \ | |
# /300\ or | 400 |
# *-----* *-----*
# ____________ ______________
boundary_edge = None
elif face_type == '003':
# -----------
# *
# / \
# / \
# *-----*
boundary_edge = None
face = list(face)
face.append(face[0])
crown_faces.append(face)
elif face_type == '004':
# ---------------
# *-----*
# | |
# | |
# *-----*
boundary_edge = None
crown_faces.append(list(face))
elif face_type == '030' or face_type == '040':
# Face is totally on the plane --> rare case...
boundary_edge = None
else:
try:
from . import mmio
mmio.write_VTP('full_debug.vtp', self.source_mesh.vertices, self.source_mesh.faces)
# mmio.write_VTP('clipped_crown_debug.vtp', clipped_crown_mesh.vertices, clipped_crown_mesh.faces)
mmio.write_VTP('crown_debug.vtp', crown_mesh.vertices, crown_mesh.faces)
except:
pass
raise Exception("Face %u clipping case %s not known." % (face_id, face_type))
# Building boundary connectivity
if boundary_edge is not None:
direct_boundary_edges[boundary_edge[0]] = boundary_edge[1]
inv_boundary_edges[boundary_edge[1]] = boundary_edge[0]
if len(intersections) > 0:
vertices = np.concatenate((vertices, intersections))
clipped_crown_mesh = Mesh(vertices, crown_faces)
# TODO: faire un merge uniquement sur la liste instersections et non sur tout le maillage clipped_crown
# FIXME: potentiellement, un bug a ete introduit ici !!! --> l'update n'est plus bon sur les dictionnaires...
new_id = clipped_crown_mesh.merge_duplicates(return_index=True, atol=1e-5) # Warning: choosing a lower value
# Updating dictionaries
direct_boundary_edges = dict(
list(zip(new_id[list(direct_boundary_edges.keys())], new_id[list(direct_boundary_edges.values())])))
inv_boundary_edges = dict(list(zip(new_id[list(inv_boundary_edges.keys())], new_id[list(inv_boundary_edges.values())])))
# Ordering boundary edges in continuous lines
closed_polygons = list()
open_lines = list()
while True:
try:
line = list()
v0_init, v1 = direct_boundary_edges.popitem()
line.append(v0_init)
line.append(v1)
v0 = v1
while True:
try:
v1 = direct_boundary_edges.pop(v0)
line.append(v1)
v0 = v1
except KeyError:
if line[0] != line[-1]:
# Trying to find an other queue
queue = list()
v0 = v0_init
while True:
try:
v1 = inv_boundary_edges[v0]
direct_boundary_edges.pop(v1)
queue.append(v1)
v0 = v1
except:
queue.reverse()
line = queue + line
# Trying to see if both end of line are not connected
pstart = clipped_crown_mesh.vertices[line[0]]
pend = clipped_crown_mesh.vertices[line[-1]]
d = np.linalg.norm(pstart-pend)
open_lines.append(line)
break
else:
closed_polygons.append(line)
break
except: # FIXME: specifier quelle exception est attendue ici !
# TODO: retirer les deux lines suivantes
if self._verbose:
print(("%u closed polygon\n%u open curve" % (len(closed_polygons), len(open_lines))))
if self._assert_closed_boundaries:
if len(open_lines) > 0:
try:
from . import mmio
mmio.write_VTP('full_debug.vtp', self.source_mesh.vertices, self.source_mesh.faces)
mmio.write_VTP('clipped_crown_debug.vtp', clipped_crown_mesh.vertices, clipped_crown_mesh.faces)
mmio.write_VTP('crown_debug.vtp', crown_mesh.vertices, crown_mesh.faces)
except:
pass
for line in open_lines:
print(line)
raise RuntimeError('Open intersection curve found with assert_closed_boundaries option enabled. Files full_debug.vtp, crown_debug.vtp and clipped_crown_debug.vtp written.')
break
output = {'clipped_crown_mesh': clipped_crown_mesh,
'closed_polygons': closed_polygons,
'open_lines': open_lines}
self.__internals__.update(output)
def _clip(self):
"""Performs clipping and assemble the clipped mesh."""
self._clip_crown_by_plane()
clipped_mesh = self.lower_mesh + self.clipped_crown_mesh
clipped_mesh.name = '_'.join((self._source_mesh.name, 'clipped'))
self.__internals__['clipped_mesh'] = clipped_mesh
return