# coding: utf-8
# /*##########################################################################
#
# Copyright (c) 2015-2017 European Synchrotron Radiation Facility
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ###########################################################################*/
"""
This module provides functions to generate indices, to check intersection
and to handle planes.
"""
from __future__ import absolute_import, division, unicode_literals
__authors__ = ["T. Vincent"]
__license__ = "MIT"
__date__ = "25/07/2016"
import logging
import numpy
from . import event
_logger = logging.getLogger(__name__)
# numpy #######################################################################
def _uniqueAlongLastAxis(a):
"""Numpy unique on the last axis of a 2D array
Implemented here as not in numpy as of writing.
See adding axis parameter to numpy.unique:
https://github.com/numpy/numpy/pull/3584/files#r6225452
:param array_like a: Input array.
:return: Unique elements along the last axis.
:rtype: numpy.ndarray
"""
assert len(a.shape) == 2
# Construct a type over last array dimension to run unique on a 1D array
if a.dtype.char in numpy.typecodes['AllInteger']:
# Bit-wise comparison of the 2 indices of a line at once
# Expect a C contiguous array of shape N, 2
uniquedt = numpy.dtype((numpy.void, a.itemsize * a.shape[-1]))
elif a.dtype.char in numpy.typecodes['Float']:
uniquedt = [('f{i}'.format(i=i), a.dtype) for i in range(a.shape[-1])]
else:
raise TypeError("Unsupported type {dtype}".format(dtype=a.dtype))
uniquearray = numpy.unique(numpy.ascontiguousarray(a).view(uniquedt))
return uniquearray.view(a.dtype).reshape((-1, a.shape[-1]))
# conversions #################################################################
[docs]def triangleToLineIndices(triangleIndices, unicity=False):
"""Generates lines indices from triangle indices.
This is generating lines indices for the edges of the triangles.
:param triangleIndices: The indices to draw a set of vertices as triangles.
:type triangleIndices: numpy.ndarray
:param bool unicity: If True remove duplicated lines,
else (the default) returns all lines.
:return: The indices to draw the edges of the triangles as lines.
:rtype: 1D numpy.ndarray of uint16 or uint32.
"""
# Makes sure indices ar packed by triangle
triangleIndices = triangleIndices.reshape(-1, 3)
# Pack line indices by triangle and by edge
lineindices = numpy.empty((len(triangleIndices), 3, 2),
dtype=triangleIndices.dtype)
lineindices[:, 0] = triangleIndices[:, :2] # edge = t0, t1
lineindices[:, 1] = triangleIndices[:, 1:] # edge =t1, t2
lineindices[:, 2] = triangleIndices[:, ::2] # edge = t0, t2
if unicity:
lineindices = _uniqueAlongLastAxis(lineindices.reshape(-1, 2))
# Make sure it is 1D
lineindices.shape = -1
return lineindices
[docs]def verticesNormalsToLines(vertices, normals, scale=1.):
"""Return vertices of lines representing normals at given positions.
:param vertices: Positions of the points.
:type vertices: numpy.ndarray with shape: (nbPoints, 3)
:param normals: Corresponding normals at the points.
:type normals: numpy.ndarray with shape: (nbPoints, 3)
:param float scale: The scale factor to apply to normals.
:returns: Array of vertices to draw corresponding lines.
:rtype: numpy.ndarray with shape: (nbPoints * 2, 3)
"""
linevertices = numpy.empty((len(vertices) * 2, 3), dtype=vertices.dtype)
linevertices[0::2] = vertices
linevertices[1::2] = vertices + scale * normals
return linevertices
[docs]def unindexArrays(mode, indices, *arrays):
"""Convert indexed GL primitives to unindexed ones.
Given indices in arrays and the OpenGL primitive they represent,
return the unindexed equivalent.
:param str mode:
Kind of primitive represented by indices.
In: points, lines, line_strip, loop, triangles, triangle_strip, fan.
:param indices: Indices in other arrays
:type indices: numpy.ndarray of dimension 1.
:param arrays: Remaining arguments are arrays to convert
:return: Converted arrays
:rtype: tuple of numpy.ndarray
"""
indices = numpy.array(indices, copy=False)
assert mode in ('points',
'lines', 'line_strip', 'loop',
'triangles', 'triangle_strip', 'fan')
if mode in ('lines', 'line_strip', 'loop'):
assert len(indices) >= 2
elif mode in ('triangles', 'triangle_strip', 'fan'):
assert len(indices) >= 3
assert indices.min() >= 0
max_index = indices.max()
for data in arrays:
assert len(data) >= max_index
if mode == 'line_strip':
unpacked = numpy.empty((2 * (len(indices) - 1),), dtype=indices.dtype)
unpacked[0::2] = indices[:-1]
unpacked[1::2] = indices[1:]
indices = unpacked
elif mode == 'loop':
unpacked = numpy.empty((2 * len(indices),), dtype=indices.dtype)
unpacked[0::2] = indices
unpacked[1:-1:2] = indices[1:]
unpacked[-1] = indices[0]
indices = unpacked
elif mode == 'triangle_strip':
unpacked = numpy.empty((3 * (len(indices) - 2),), dtype=indices.dtype)
unpacked[0::3] = indices[:-2]
unpacked[1::3] = indices[1:-1]
unpacked[2::3] = indices[2:]
indices = unpacked
elif mode == 'fan':
unpacked = numpy.empty((3 * (len(indices) - 2),), dtype=indices.dtype)
unpacked[0::3] = indices[0]
unpacked[1::3] = indices[1:-1]
unpacked[2::3] = indices[2:]
indices = unpacked
return tuple(numpy.ascontiguousarray(data[indices]) for data in arrays)
[docs]def trianglesNormal(positions):
"""Return normal for each triangle.
:param positions: Serie of triangle's corners
:type positions: numpy.ndarray of shape (NbTriangles*3, 3)
:return: Normals corresponding to each position.
:rtype: numpy.ndarray of shape (NbTriangles, 3)
"""
assert positions.ndim == 2
assert positions.shape[1] == 3
positions = numpy.array(positions, copy=False).reshape(-1, 3, 3)
normals = numpy.cross(positions[:, 1] - positions[:, 0],
positions[:, 2] - positions[:, 0])
# Normalize normals
if numpy.version.version < '1.8.0':
# Debian 7 support: numpy.linalg.norm has no axis argument
norms = numpy.array(tuple(numpy.linalg.norm(vec) for vec in normals),
dtype=normals.dtype)
else:
norms = numpy.linalg.norm(normals, axis=1)
norms[norms == 0] = 1
return normals / norms.reshape(-1, 1)
# grid ########################################################################
[docs]def gridVertices(dim0Array, dim1Array, dtype):
"""Generate an array of 2D positions from 2 arrays of 1D coordinates.
:param dim0Array: 1D array-like of coordinates along the first dimension.
:param dim1Array: 1D array-like of coordinates along the second dimension.
:param numpy.dtype dtype: Data type of the output array.
:return: Array of grid coordinates.
:rtype: numpy.ndarray with shape: (len(dim0Array), len(dim1Array), 2)
"""
grid = numpy.empty((len(dim0Array), len(dim1Array), 2), dtype=dtype)
grid.T[0, :, :] = dim0Array
grid.T[1, :, :] = numpy.array(dim1Array, copy=False)[:, None]
return grid
[docs]def triangleStripGridIndices(dim0, dim1):
"""Generate indices to draw a grid of vertices as a triangle strip.
Vertices are expected to be stored as row-major (i.e., C contiguous).
:param int dim0: The number of rows of vertices.
:param int dim1: The number of columns of vertices.
:return: The vertex indices
:rtype: 1D numpy.ndarray of uint32
"""
assert dim0 >= 2
assert dim1 >= 2
# Filling a row of squares +
# an index before and one after for degenerated triangles
indices = numpy.empty((dim0 - 1, 2 * (dim1 + 1)), dtype=numpy.uint32)
# Init indices with minimum indices for each row of squares
indices[:] = (dim1 * numpy.arange(dim0 - 1, dtype=numpy.uint32))[:, None]
# Update indices with offset per row of squares
offset = numpy.arange(dim1, dtype=numpy.uint32)
indices[:, 1:-1:2] += offset
offset += dim1
indices[:, 2::2] += offset
indices[:, -1] += offset[-1]
# Remove extra indices for degenerated triangles before returning
return indices.ravel()[1:-1]
# Alternative:
# indices = numpy.zeros(2 * dim1 * (dim0 - 1) + 2 * (dim0 - 2),
# dtype=numpy.uint32)
#
# offset = numpy.arange(dim1, dtype=numpy.uint32)
# for d0Index in range(dim0 - 1):
# start = 2 * d0Index * (dim1 + 1)
# end = start + 2 * dim1
# if d0Index != 0:
# indices[start - 2] = offset[-1]
# indices[start - 1] = offset[0]
# indices[start:end:2] = offset
# offset += dim1
# indices[start + 1:end:2] = offset
# return indices
[docs]def linesGridIndices(dim0, dim1):
"""Generate indices to draw a grid of vertices as lines.
Vertices are expected to be stored as row-major (i.e., C contiguous).
:param int dim0: The number of rows of vertices.
:param int dim1: The number of columns of vertices.
:return: The vertex indices.
:rtype: 1D numpy.ndarray of uint32
"""
# Horizontal and vertical lines
nbsegmentalongdim1 = 2 * (dim1 - 1)
nbsegmentalongdim0 = 2 * (dim0 - 1)
indices = numpy.empty(nbsegmentalongdim1 * dim0 +
nbsegmentalongdim0 * dim1,
dtype=numpy.uint32)
# Line indices over dim0
onedim1line = (numpy.arange(nbsegmentalongdim1,
dtype=numpy.uint32) + 1) // 2
indices[:dim0 * nbsegmentalongdim1] = \
(dim1 * numpy.arange(dim0, dtype=numpy.uint32)[:, None] +
onedim1line[None, :]).ravel()
# Line indices over dim1
onedim0line = (numpy.arange(nbsegmentalongdim0,
dtype=numpy.uint32) + 1) // 2
indices[dim0 * nbsegmentalongdim1:] = \
(numpy.arange(dim1, dtype=numpy.uint32)[:, None] +
dim1 * onedim0line[None, :]).ravel()
return indices
# intersection ################################################################
[docs]def angleBetweenVectors(refVector, vectors, norm=None):
"""Return the angle between 2 vectors.
:param refVector: Coordinates of the reference vector.
:type refVector: numpy.ndarray of shape: (NCoords,)
:param vectors: Coordinates of the vector(s) to get angle from reference.
:type vectors: numpy.ndarray of shape: (NCoords,) or (NbVector, NCoords)
:param norm: A direction vector giving an orientation to the angles
or None.
:returns: The angles in radians in [0, pi] if norm is None
else in [0, 2pi].
:rtype: float or numpy.ndarray of shape (NbVectors,)
"""
singlevector = len(vectors.shape) == 1
if singlevector: # Make it a 2D array for the computation
vectors = vectors.reshape(1, -1)
assert len(refVector.shape) == 1
assert len(vectors.shape) == 2
assert len(refVector) == vectors.shape[1]
# Normalize vectors
refVector /= numpy.linalg.norm(refVector)
vectors = numpy.array([v / numpy.linalg.norm(v) for v in vectors])
dots = numpy.sum(refVector * vectors, axis=-1)
angles = numpy.arccos(numpy.clip(dots, -1., 1.))
if norm is not None:
signs = numpy.sum(norm * numpy.cross(refVector, vectors), axis=-1) < 0.
angles[signs] = numpy.pi * 2. - angles[signs]
return angles[0] if singlevector else angles
[docs]def segmentPlaneIntersect(s0, s1, planeNorm, planePt):
"""Compute the intersection of a segment with a plane.
:param s0: First end of the segment
:type s0: 1D numpy.ndarray-like of length 3
:param s1: Second end of the segment
:type s1: 1D numpy.ndarray-like of length 3
:param planeNorm: Normal vector of the plane.
:type planeNorm: numpy.ndarray of shape: (3,)
:param planePt: A point of the plane.
:type planePt: numpy.ndarray of shape: (3,)
:return: The intersection points. The number of points goes
from 0 (no intersection) to 2 (segment in the plane)
:rtype: list of numpy.ndarray
"""
s0, s1 = numpy.asarray(s0), numpy.asarray(s1)
segdir = s1 - s0
dotnormseg = numpy.dot(planeNorm, segdir)
if dotnormseg == 0:
# line and plane are parallels
if numpy.dot(planeNorm, planePt - s0) == 0: # segment is in plane
return [s0, s1]
else: # No intersection
return []
alpha = - numpy.dot(planeNorm, s0 - planePt) / dotnormseg
if 0. <= alpha <= 1.: # Intersection with segment
return [s0 + alpha * segdir]
else: # intersection outside segment
return []
[docs]def boxPlaneIntersect(boxVertices, boxLineIndices, planeNorm, planePt):
"""Return intersection points between a box and a plane.
:param boxVertices: Position of the corners of the box.
:type boxVertices: numpy.ndarray with shape: (8, 3)
:param boxLineIndices: Indices of the box edges.
:type boxLineIndices: numpy.ndarray-like with shape: (12, 2)
:param planeNorm: Normal vector of the plane.
:type planeNorm: numpy.ndarray of shape: (3,)
:param planePt: A point of the plane.
:type planePt: numpy.ndarray of shape: (3,)
:return: The found intersection points
:rtype: numpy.ndarray with 2 dimensions
"""
segments = numpy.take(boxVertices, boxLineIndices, axis=0)
points = set() # Gather unique intersection points
for seg in segments:
for point in segmentPlaneIntersect(seg[0], seg[1], planeNorm, planePt):
points.add(tuple(point))
points = numpy.array(list(points))
if len(points) <= 2:
return numpy.array(())
elif len(points) == 3:
return points
else: # len(points) > 3
# Order point to have a polyline lying on the unit cube's faces
vectors = points - numpy.mean(points, axis=0)
angles = angleBetweenVectors(vectors[0], vectors, planeNorm)
points = numpy.take(points, numpy.argsort(angles), axis=0)
return points
# Plane #######################################################################
[docs]class Plane(event.Notifier):
"""Object handling a plane and notifying plane changes.
:param point: A point on the plane.
:type point: 3-tuple of float.
:param normal: Normal of the plane.
:type normal: 3-tuple of float.
"""
def __init__(self, point=(0., 0., 0.), normal=(0., 0., 1.)):
super(Plane, self).__init__()
assert len(point) == 3
self._point = numpy.array(point, copy=True, dtype=numpy.float32)
assert len(normal) == 3
self._normal = numpy.array(normal, copy=True, dtype=numpy.float32)
self.notify()
[docs] def setPlane(self, point=None, normal=None):
"""Set plane point and normal and notify.
:param point: A point on the plane.
:type point: 3-tuple of float or None.
:param normal: Normal of the plane.
:type normal: 3-tuple of float or None.
"""
planechanged = False
if point is not None:
assert len(point) == 3
point = numpy.array(point, copy=True, dtype=numpy.float32)
if not numpy.all(numpy.equal(self._point, point)):
self._point = point
planechanged = True
if normal is not None:
assert len(normal) == 3
normal = numpy.array(normal, copy=True, dtype=numpy.float32)
norm = numpy.linalg.norm(normal)
if norm != 0.:
normal /= norm
if not numpy.all(numpy.equal(self._normal, normal)):
self._normal = normal
planechanged = True
if planechanged:
_logger.debug('Plane updated:\n\tpoint: %s\n\tnormal: %s',
str(self._point), str(self._normal))
self.notify()
@property
def point(self):
"""A point on the plane."""
return self._point.copy()
@point.setter
def point(self, point):
self.setPlane(point=point)
@property
def normal(self):
"""The (normalized) normal of the plane."""
return self._normal.copy()
@normal.setter
def normal(self, normal):
self.setPlane(normal=normal)
@property
def parameters(self):
"""Plane equation parameters: a*x + b*y + c*z + d = 0."""
return numpy.append(self._normal,
- numpy.dot(self._point, self._normal))
@parameters.setter
def parameters(self, parameters):
assert len(parameters) == 4
parameters = numpy.array(parameters, dtype=numpy.float32)
# Normalize normal
norm = numpy.linalg.norm(parameters[:3])
if norm != 0:
parameters /= norm
normal = parameters[:3]
point = - parameters[3] * normal
self.setPlane(point, normal)
@property
def isPlane(self):
"""True if a plane is defined (i.e., ||normal|| != 0)."""
return numpy.any(self.normal != 0.)
[docs] def move(self, step):
"""Move the plane of step along the normal."""
self.point += step * self.normal