# coding: utf-8
# /*##########################################################################
#
# Copyright (c) 2015-2020 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 triangleStripToTriangles(strip):
"""Convert a triangle strip to a set of triangles.
The order of the corners is inverted for odd triangles.
:param numpy.ndarray strip:
Array of triangle corners of shape (N, 3).
N must be at least 3.
:return: Equivalent triangles corner as an array of shape (N, 3, 3)
:rtype: numpy.ndarray
"""
strip = numpy.array(strip).reshape(-1, 3)
assert len(strip) >= 3
triangles = numpy.empty((len(strip) - 2, 3, 3), dtype=strip.dtype)
triangles[0::2, 0] = strip[0:-2:2]
triangles[0::2, 1] = strip[1:-1:2]
triangles[0::2, 2] = strip[2::2]
triangles[1::2, 0] = strip[3::2]
triangles[1::2, 1] = strip[2:-1:2]
triangles[1::2, 2] = strip[1:-2:2]
return triangles
[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
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
[docs]def clipSegmentToBounds(segment, bounds):
"""Clip segment to volume aligned with axes.
:param numpy.ndarray segment: (p0, p1)
:param numpy.ndarray bounds: (lower corner, upper corner)
:return: Either clipped (p0, p1) or None if outside volume
:rtype: Union[None,List[numpy.ndarray]]
"""
segment = numpy.array(segment, copy=False)
bounds = numpy.array(bounds, copy=False)
p0, p1 = segment
# Get intersection points of ray with volume boundary planes
# Line equation: P = offset * delta + p0
delta = p1 - p0
deltaNotZero = numpy.array(delta, copy=True)
deltaNotZero[deltaNotZero == 0] = numpy.nan # Invalidated to avoid division by zero
offsets = ((bounds - p0) / deltaNotZero).reshape(-1)
points = offsets.reshape(-1, 1) * delta + p0
# Avoid precision errors by using bounds value
points.shape = 2, 3, 3 # Reshape 1 point per bound value
for dim in range(3):
points[:, dim, dim] = bounds[:, dim]
points.shape = -1, 3 # Set back to 2D array
# Find intersection points that are included in the volume
mask = numpy.logical_and(numpy.all(bounds[0] <= points, axis=1),
numpy.all(points <= bounds[1], axis=1))
intersections = numpy.unique(offsets[mask])
if len(intersections) != 2:
return None
intersections.sort()
# Do p1 first as p0 is need to compute it
if intersections[1] < 1: # clip p1
segment[1] = intersections[1] * delta + p0
if intersections[0] > 0: # clip p0
segment[0] = intersections[0] * delta + p0
return segment
[docs]def segmentVolumeIntersect(segment, nbins):
"""Get bin indices intersecting with segment
It should work with N dimensions.
Coordinate convention (z, y, x) or (x, y, z) should not matter
as long as segment and nbins are consistent.
:param numpy.ndarray segment:
Segment end points as a 2xN array of coordinates
:param numpy.ndarray nbins:
Shape of the volume with same coordinates order as segment
:return: List of bins indices as a 2D array or None if no bins
:rtype: Union[None,numpy.ndarray]
"""
segment = numpy.asarray(segment)
nbins = numpy.asarray(nbins)
assert segment.ndim == 2
assert segment.shape[0] == 2
assert nbins.ndim == 1
assert segment.shape[1] == nbins.size
dim = len(nbins)
bounds = numpy.array((numpy.zeros_like(nbins), nbins))
segment = clipSegmentToBounds(segment, bounds)
if segment is None:
return None # Segment outside volume
p0, p1 = segment
# Get intersections
# Get coordinates of bin edges crossing the segment
clipped = numpy.ceil(numpy.clip(segment, 0, nbins))
start = numpy.min(clipped, axis=0)
stop = numpy.max(clipped, axis=0) # stop is NOT included
edgesByDim = [numpy.arange(start[i], stop[i]) for i in range(dim)]
# Line equation: P = t * delta + p0
delta = p1 - p0
# Get bin edge/line intersections as sorted points along the line
# Get corresponding line parameters
t = []
if numpy.all(0 <= p0) and numpy.all(p0 <= nbins):
t.append([0.]) # p0 within volume, add it
t += [(edgesByDim[i] - p0[i]) / delta[i] for i in range(dim) if delta[i] != 0]
if numpy.all(0 <= p1) and numpy.all(p1 <= nbins):
t.append([1.]) # p1 within volume, add it
t = numpy.concatenate(t)
t.sort(kind='mergesort')
# Remove duplicates
unique = numpy.ones((len(t),), dtype=bool)
numpy.not_equal(t[1:], t[:-1], out=unique[1:])
t = t[unique]
if len(t) < 2:
return None # Not enough intersection points
# bin edges/line intersection points
points = t.reshape(-1, 1) * delta + p0
centers = (points[:-1] + points[1:]) / 2.
bins = numpy.floor(centers).astype(numpy.int64)
return bins
# 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
[docs] def segmentIntersection(self, s0, s1):
"""Compute the plane intersection with segment [s0, s1].
: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
:return: The intersection points. The number of points goes
from 0 (no intersection) to 2 (segment in the plane)
:rtype: list of 1D numpy.ndarray
"""
if not self.isPlane:
return []
else:
return segmentPlaneIntersect(s0, s1, self.normal, self.point)
