Commit 5069715e authored by Daniel Boeckenhoff's avatar Daniel Boeckenhoff
Browse files

tfields.Mesh3D cutting works - except for templates on fields

parent 07604102
......@@ -388,6 +388,9 @@ class Tensors(AbstractNdarray):
tensors = np.empty(tensors.shape, dtype=tensors.dtype)
elif dim is not None:
tensors = np.empty((0, dim))
if issubclass(type(tensors), np.ndarray):
# np.empty
pass
else:
raise ValueError("Empty tensors need dimension "
"parameter 'dim'.")
......@@ -1454,10 +1457,11 @@ class TensorMaps(TensorFields):
dims.append(map_field.dim)
else:
maps[mp_idx] = TensorFields.merged(maps[mp_idx], map_field)
kwargs['maps'] = maps
# kwargs['maps'] = maps
obj = super(TensorMaps, cls).merged(*objects, **kwargs)
return cls.__new__(cls, obj, maps=maps)
inst = super(TensorMaps, cls).merged(*objects, **kwargs)
inst = cls.__new__(cls, inst, maps=maps)
return inst
def equal(self, other, **kwargs):
"""
......
......@@ -108,3 +108,4 @@ else:
from .grid import igrid
from . import stats
from .stats import mode, median, mean
from . import symbolics
......@@ -15,7 +15,6 @@ import ioTools
import mplTools
import decoTools
import pyTools
import symTools
from sympy.abc import y, z
from scipy.spatial import ConvexHull
import matplotlib.pyplot as plt
......@@ -26,6 +25,118 @@ import cuttingTree
logger = loggingTools.Logger(__name__)
def _dist_from_plane(point, plane):
return plane['normal'].dot(point) + plane['d']
def _segment_plane_intersection(p0, p1, plane):
"""
Returns:
points, direction
"""
distance0 = _dist_from_plane(p0, plane)
distance1 = _dist_from_plane(p1, plane)
p0OnPlane = abs(distance0) < np.finfo(float).eps
p1OnPlane = abs(distance1) < np.finfo(float).eps
outPoints = []
direction = 0
if p0OnPlane:
outPoints.append(p0)
if p1OnPlane:
outPoints.append(p1)
# remove duplicate points
if len(outPoints) > 1:
outPoints = tfields.Points3D(np.unique(outPoints, axis=0))
if p0OnPlane and p1OnPlane:
return outPoints, direction
if distance0 * distance1 > np.finfo(float).eps:
return outPoints, direction
direction = np.sign(distance0)
if abs(distance0) < np.finfo(float).eps:
return outPoints, direction
elif abs(distance1) < np.finfo(float).eps:
return outPoints, direction
if abs(distance0 - distance1) > np.finfo(float).eps:
t = distance0 / (distance0 - distance1)
else:
return outPoints, direction
outPoints.append(p0 + t * (p1 - p0))
# remove duplicate points
if len(outPoints) > 1:
outPoints = tfields.Points3D(np.unique(outPoints, axis=0))
return outPoints, direction
def _intersect(triangle, plane, vertices_rejected):
"""
Intersect a triangle with a plane. Give the info, which side of the
triangle is rejected by passing the mask vertices_rejected
Returns:
tuple: points, faces
"""
nTrue = vertices_rejected.count(True)
lonely_bool = True if nTrue == 1 else False
index = vertices_rejected.index(lonely_bool)
s0, d0 = _segment_plane_intersection(triangle[0], triangle[1], plane)
s1, d1 = _segment_plane_intersection(triangle[1], triangle[2], plane)
s2, d2 = _segment_plane_intersection(triangle[2], triangle[0], plane)
singlePoint = triangle[index]
couplePoints = [triangle[j] for j in range(3)
if not vertices_rejected[j] == lonely_bool]
# TODO handle special cases. For now triangles with at least two points on plane are excluded
new_points = None
faces = None
if len(s0) == 2:
# both points on plane
return None, None
if len(s1) == 2:
# both points on plane
return None, None
if len(s2) == 2:
# both points on plane
return None, None
if lonely_bool:
# one new triangle
if len(s0) == 1 and len(s1) == 1:
new_points = np.array([couplePoints[0], couplePoints[1], s0[0], s1[0]])
faces = np.array([[0, 2, 1], [1, 2, 3]])
elif len(s1) == 1 and len(s2) == 1:
new_points = np.array([couplePoints[0], couplePoints[1], s1[0], s2[0]])
faces = np.array([[0, 1, 2], [0, 2, 3]])
elif len(s0) == 1 and len(s2) == 1:
new_points = np.array([couplePoints[0], couplePoints[1], s0[0], s2[0]])
faces = np.array([[0, 1, 2], [1, 3, 2]])
else:
return None, None
else:
# two new triangles
if len(s0) == 1 and len(s1) == 1:
new_points = np.array([singlePoint, s0[0], s1[0]])
faces = np.array([[0, 2, 1]])
elif len(s1) == 1 and len(s2) == 1:
new_points = np.array([singlePoint, s1[0], s2[0]])
faces = np.array([[0, 2, 1]])
elif len(s0) == 1 and len(s2) == 1:
new_points = np.array([singlePoint, s0[0], s2[0]])
faces = np.array([[0, 1, 2]])
else:
return None, None
return new_points, faces
def scalars_to_fields(scalars):
scalars = np.array(scalars)
if len(scalars.shape) == 1:
......@@ -109,10 +220,9 @@ class Mesh3D(tfields.TensorMaps):
array([[0, 1, 2]])
"""
def __new__(cls, tensors, **kwargs):
def __new__(cls, tensors, *fields, **kwargs):
if not issubclass(type(tensors), Mesh3D):
kwargs['dim'] = 3
fields = kwargs.pop('fields', [])
faces = kwargs.pop('faces', None)
faceScalars = kwargs.pop('faceScalars', [])
maps = kwargs.pop('maps', None)
......@@ -647,156 +757,18 @@ class Mesh3D(tfields.TensorMaps):
True
"""
# log = logger.new()
faceIndices = np.arange(self.faces.shape[0])
cents = tfields.Points3D(sub_mesh.getCentroids())
scalars = []
# nCents = cents.shape[0]
# for i in range(nCents):
# faceMask = self.pointsInMesh(cents[i: i + 1], delta=delta)[0]
# if True not in faceMask:
# log.warning("No point face was assigned to this. I will retry.")
# faceMask = self.pointsInMesh(cents[i: i + 1], delta=delta * 100)[0]
# if True not in faceMask:
# raise ValueError()
# else:
# log.info("Centroid {i} / {nCents}".format(**locals()))
# scalars.append(faceIndices[faceMask])
mask = self.pointsInMesh(cents, delta=delta)
scalars = [faceIndices[faceMask] for faceMask in mask]
inst = sub_mesh.copy()
inst.setScalarArray(0, scalars)
return inst
def _distFromPlane(self, point, plane):
return plane['normal'].dot(point) + plane['d']
def _getSegmentPlaneIntersection(self, p0, p1, plane):
distance0 = self._distFromPlane(p0, plane)
distance1 = self._distFromPlane(p1, plane)
p0OnPlane = abs(distance0) < np.finfo(float).eps
p1OnPlane = abs(distance1) < np.finfo(float).eps
outPoints = []
direction = 0
if p0OnPlane:
outPoints.append(p0)
if p1OnPlane:
outPoints.append(p1)
# remove duplicate points
if len(outPoints) > 1:
outPoints = tfields.Points3D(np.unique(outPoints, axis=0))
if p0OnPlane and p1OnPlane:
return outPoints, direction
if distance0 * distance1 > np.finfo(float).eps:
return outPoints, direction
direction = np.sign(distance0)
if abs(distance0) < np.finfo(float).eps:
return outPoints, direction
elif abs(distance1) < np.finfo(float).eps:
return outPoints, direction
if abs(distance0 - distance1) > np.finfo(float).eps:
t = distance0 / (distance0 - distance1)
else:
return outPoints, direction
outPoints.append(p0 + t * (p1 - p0))
# remove duplicate points
if len(outPoints) > 1:
outPoints = tfields.Points3D(np.unique(outPoints, axis=0))
return outPoints, direction
def _intersect(self, triangle, plane, facePointsRejected):
nTrue = facePointsRejected.count(True)
lonely_bool = True if nTrue == 1 else False
index = facePointsRejected.index(lonely_bool)
s0, d0 = self._getSegmentPlaneIntersection(triangle[0], triangle[1], plane)
s1, d1 = self._getSegmentPlaneIntersection(triangle[1], triangle[2], plane)
s2, d2 = self._getSegmentPlaneIntersection(triangle[2], triangle[0], plane)
singlePoint = triangle[index]
couplePoints = [triangle[j] for j in range(3)
if not facePointsRejected[j] == lonely_bool]
# TODO handle special cases. For now triangles with at least two points on plane are excluded
new_points = None
faces = None
if len(s0) == 2:
# both points on plane
return None, None
if len(s1) == 2:
# both points on plane
return None, None
if len(s2) == 2:
# both points on plane
return None, None
if lonely_bool:
# one new triangle
if len(s0) == 1 and len(s1) == 1:
new_points = np.array([couplePoints[0], couplePoints[1], s0[0], s1[0]])
faces = np.array([[0, 2, 1], [1, 2, 3]])
elif len(s1) == 1 and len(s2) == 1:
new_points = np.array([couplePoints[0], couplePoints[1], s1[0], s2[0]])
faces = np.array([[0, 1, 2], [0, 2, 3]])
elif len(s0) == 1 and len(s2) == 1:
new_points = np.array([couplePoints[0], couplePoints[1], s0[0], s2[0]])
faces = np.array([[0, 1, 2], [1, 3, 2]])
else:
return None, None
else:
# two new triangles
if len(s0) == 1 and len(s1) == 1:
new_points = np.array([singlePoint, s0[0], s1[0]])
faces = np.array([[0, 2, 1]])
elif len(s1) == 1 and len(s2) == 1:
new_points = np.array([singlePoint, s1[0], s2[0]])
faces = np.array([[0, 2, 1]])
elif len(s0) == 1 and len(s2) == 1:
new_points = np.array([singlePoint, s0[0], s2[0]])
faces = np.array([[0, 1, 2]])
else:
return None, None
return new_points, faces
def _cut_sympy(self, expression, at_intersection="remove", _template=None):
def _cut_sympy(self, expression, at_intersection="remove", _in_recursion=False):
"""
Partition the mesh with the cuts given
Args:
expression (sympy logical expression): see Tensors.cut
at_intersection (str): option switch for expression (see
tfields.Mesh3D.cut)
Examples:
Cut away the first face / select the last face only
>>> m = tfields.Mesh3D([[1,2,3], [3,3,3], [0,0,0], [5,6,7]],
... faces=[[0, 1, 2], [1, 2, 3]],
... faceScalars=[[1,2], [6,7]])
>>> mNew = tfields.Mesh3D([[3,3,3], [0,0,0], [5,6,7]], [[0, 1, 2]], faceScalars=[[1]])
>>> mCut = m.cut(mNew)
>>> mCut.faceScalars
array([[ 6., 7.]])
Cut with 'return_template=True' will return an instruction
on how to conduct the cut fast (template)
>>> from sympy.abc import y
>>> m = tfields.Mesh3D([[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,2,0], [1,2,0]],
... faces=[[0,1,2],[2,3,0],[3,2,5],[5,4,3]],
... faceScalars=[[1],[2],[3],[4]])
>>> mCut, template = m.cut(y < 1.5, at_intersection='split',
... return_template=True)
>>> mCut.faceScalars.T
array([[ 1., 2., 3., 3., 4.]])
Applying the template results in the same result as applying the cut
with the condition
>>> mCut2 = m.cut(template)
>>> bool((mCut == mCut2).all())
True
>>> bool((mCut.faceScalars == mCut2.faceScalars).all())
True
Partition the mesh with the cuts given and return the template
"""
eps = 0.000000001
......@@ -804,27 +776,35 @@ class Mesh3D(tfields.TensorMaps):
if len(self) == 0:
return self.copy(), self.copy()
inst = self.copy()
'''
add the indices of the vertices and maps to the fields. They will be
removed afterwards
'''
if not _in_recursion:
inst.fields.append(tfields.Tensors(np.arange(len(inst))))
for mp in inst.maps:
mp.fields.append(tfields.Tensors(np.arange(len(mp))))
# mask for points that do not fulfill the cut expression
mask = self.evalf(expression)
mask = inst.evalf(expression)
# remove the points
inst = self[mask].copy()
faceKeepMask = self.getFaceMask(mask)
scalarMap = np.arange(self.nfaces())[faceKeepMask]
scalarMap = scalarMap.reshape((-1, 1)).astype(float)
if not any(~mask):
# no faces have been removed.
return inst, Mesh3D(inst, faceScalars=scalarMap)
if all(~mask):
# all faces have been removed)
return inst, inst.copy()
# add points and faces intersecting with the plane
if at_intersection == 'split' or at_intersection == 'splitRough':
expression_parts = symTools.getExpressionParts(expression) # TODO: move to tfields lib
inst_pre_mask = inst.copy()
inst = inst[mask]
if not any(~mask) or all(~mask):
pass
elif at_intersection == 'split' or at_intersection == 'splitRough':
# add points and faces intersecting with the plane
expression_parts = tfields.lib.symbolics.split_expression(expression)
'''
define a new mesh that will be merged with the existing one
the new mesh will describe the faces that are at the border of the
cuts
'''
if len(expression_parts) > 1:
new_mesh = self.copy()
new_mesh = inst_pre_mask.copy()
if at_intersection == 'splitRough':
"""
the following is, to speed up the process. Problem is, that
......@@ -835,21 +815,22 @@ class Mesh3D(tfields.TensorMaps):
"""
faceIntersMask = np.full((self.faces.shape[0]), False, dtype=bool)
for i, face in enumerate(self.faces):
facePointsRejected = [-mask[f] for f in face]
face_on_edge = any(facePointsRejected) and not all(facePointsRejected)
vertices_rejected = [-mask[f] for f in face]
face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
if face_on_edge:
faceIntersMask[i] = True
new_mesh.removeFaces(-faceIntersMask)
for exprPart in expression_parts:
new_mesh, _template = new_mesh._cut_sympy(exprPart,
at_intersection='split',
_template=_template)
new_mesh, _ = new_mesh._cut_sympy(exprPart,
at_intersection='split',
_in_recursion=True)
elif len(expression_parts) == 1:
# TODO maps[0] -> smthng like self.get_map(dim=3)
points = [sympy.symbols('x0, y0, z0'),
sympy.symbols('x1, y1, z1'),
sympy.symbols('x2, y2, z2')]
plane_sympy = symTools.getPlane(expression)
plane_sympy = tfields.lib.symbolics.to_plane(expression)
norm_sympy = np.array(plane_sympy.normal_vector).astype(float)
d = -norm_sympy.dot(np.array(plane_sympy.p1).astype(float))
plane = {'normal': norm_sympy, 'd': d}
......@@ -857,24 +838,28 @@ class Mesh3D(tfields.TensorMaps):
norm_vectors = self.triangles.norms()
new_points = np.empty((0, 3))
new_faces = np.empty((0, 3))
new_map_fields = [[] for field in self.maps[0].fields]
new_fields = [tfields.Tensors(np.empty((0,) + field.shape[1:]),
coordSys=field.coordSys)
for field in inst.fields]
new_map_fields = [[] for field in inst.maps[0].fields]
new_norm_vectors = []
newScalarMap = []
n_new = 0
for i, face in enumerate(self.faces):
for i, face in enumerate(self.maps[0]):
"""
facePointsRejected is a mask for each face that is True, where
vertices_rejected is a mask for each face that is True, where
a Point is on the rejected side of the plane
"""
facePointsRejected = [~mask[f] for f in face]
face_on_edge = any(facePointsRejected) and not all(facePointsRejected)
vertices_rejected = [~mask[f] for f in face]
face_on_edge = any(vertices_rejected) and not all(vertices_rejected)
if face_on_edge:
"""
define the two lines that are intersecting the plane:
one point will be alone on one side of the cutting plane (singlePoint)
the other two will be on ther other side (couplePoints)
"""
nTrue = facePointsRejected.count(True)
nTrue = vertices_rejected.count(True)
# the bool that occures on
lonely_bool = True if nTrue == 1 else False
......@@ -882,66 +867,93 @@ class Mesh3D(tfields.TensorMaps):
"""
Add the intersection points and faces
"""
# tick = time()
if lonely_bool:
# singlePoint on cut away side
newP, newF = self._intersect(triangle_points, plane, facePointsRejected)
if newP is not None:
print "______________"
print triangle_points, vertices_rejected
newP, newF = _intersect(triangle_points, plane, vertices_rejected)
print newP, newF
if newP is not None:
if lonely_bool:
'''
singlePoint on cut away side
build two new triangles with two known and two new
points
'''
newP = np.array(newP)
np.place(newP, np.logical_and(0 < newP, newP <= eps), 0.0)
new_points = np.concatenate((new_points, newP))
new_fields = [tfields.Tensors.merged(field,
np.full((4,) + field.shape[1:], np.nan))
for field in new_fields]
new_faces = np.concatenate((new_faces, newF +
n_new))
for fld_idx, fld in enumerate(self.maps[0].fields):
new_map_fields[fld_idx].extend([fld[i]] * 2)
for field_idx, field in enumerate(self.maps[0].fields):
new_map_fields[field_idx].extend([field[i]] * 2)
if not _in_recursion:
new_map_fields[-1].extend([i] * 2)
new_norm_vectors.extend([norm_vectors[i]] * 2)
newScalarMap.extend([i] * 2)
n_new += 4
else:
pass
else:
# singlePoint on keep side
newP, newF = self._intersect(triangle_points, plane, facePointsRejected)
if newP is not None:
'''
singlePoint on keep side
build one new triangle with one known and two new points
'''
newP = np.array(newP)
np.place(newP, np.logical_and(0 < newP, newP <= eps), 0.0)
new_points = np.concatenate((new_points, newP))
new_fields = [tfields.Tensors.merged(field,
np.full((3,) + field.shape[1:], np.nan))
for field in new_fields]
new_faces = np.concatenate((new_faces, newF +
n_new))
for fld_idx, fld in enumerate(self.maps[0].fields):
new_map_fields[fld_idx].extend([fld[i]] * 1)
for field_idx, field in enumerate(self.maps[0].fields):
new_map_fields[field_idx].extend([field[i]] * 1)
if not _in_recursion:
new_map_fields[-1].extend([i] * 1)
new_norm_vectors.extend([norm_vectors[i]] * 1)
newScalarMap.extend([i])
n_new += 3
else:
pass
new_points = [[float(p[0]), float(p[1]), float(p[2])] for p in new_points] # TODO:[map(float, p) for p in new_points]
face_map = tfields.TensorFields(new_faces, *new_map_fields, dtype=int,
new_points = [map(float, p) for p in new_points]
face_map = tfields.TensorFields(new_faces, *new_map_fields,
dtype=int,
coordSys=inst.coordSys)
new_mesh = self.__class__(new_points, maps=[face_map],
coordSys=inst.coordSys,
faceScalars=new_map_fields)
new_mesh = tfields.Mesh3D(new_points,
*new_fields,
maps=[face_map],
coordSys=inst.coordSys)
new_mesh.align_norms(new_norm_vectors)
newScalarMap = np.array(newScalarMap).reshape((-1, 1)).astype(float)
scalarMap = np.concatenate((scalarMap, newScalarMap))
else:
raise ValueError("Sympy expression is not splitable.")
'''
Merge parts of mesh that are clearly in the cuts (inst) with
parts of mesh on the edge of the cuts, where new vertices and faces
parts of mesh on the edge of the cuts, where new vertices and faces
have been defined
'''
inst = self.__class__.merged(inst, new_mesh)
inst = tfields.Mesh3D.merged(inst, new_mesh)
elif at_intersection == 'remove':
pass
else:
raise AttributeError("No at_intersection method called {at_intersection} "
"implemented".format(**locals()))
template = inst.copy()
template.faceScalars = scalarMap
if _in_recursion:
template = None
else:
print inst.fields
template_field = inst.fields.pop(-1)
template_maps = []
for mp in inst.maps:
t_mp = tfields.TensorFields(tfields.Tensors(mp),
mp.fields.pop(-1))
template_maps.append(t_mp)
template = tfields.Mesh3D(tfields.Tensors(inst),
template_field,
maps=template_maps)
print template.fields
return inst, template
def _cut_template(self, template):
......@@ -953,21 +965,22 @@ class Mesh3D(tfields.TensorMaps):
>>> import tfields
>>> import numpy as np
Mesh
Build mesh
>>> mmap = tfields.TensorFields([[0, 1, 2], [0, 3, 4]],
... [[42, 21], [-42, -21]])
>>> m = tfields.Mesh3D([[0]*3, [1]*3, [2]*3, [3]*3, [4]*3],
... fields=[[0.0, 0.1, 0.2, 0.3, 0.4],
... [0.0, -0.1, -0.2, -0.3, -0.4]],
... [0.0, 0.1, 0.2, 0.3, 0.4],
... [0.0, -0.1, -0.2, -0.3, -0.4],
... maps=[mmap])
Template
Build template
>>> tmap = tfields.TensorFields([[0, 3, 4], [0, 1, 2]],
... [1, 0])
>>> t = tfields.Mesh3D([[0]*3, [-1]*3, [-2]*3, [-3]*3, [-4]*3],
... fields=[[1, 0, 3, 2, 4]],
... [1, 0, 3, 2, 4]
... maps=[tmap])
Use template as instruction to make a fast cut
>>> res = m._cut_template(t)
>>> assert np.array_equal(res.fields,
... [[0.1, 0.0, 0.3, 0.2, 0.4],
......@@ -977,11 +990,14 @@ class Mesh3D(tfields.TensorMaps):
... [[-42, -21], [42, 21]])
"""
# Redirect fields
if template.fields:
fields = [field[template.fields[0].astype(int)]
for field in self.fields]
else:
fields = []
# Redirect maps fields
maps = []
for mp, template_mp in zip(self.maps, template.maps):