diff options
Diffstat (limited to 'src/silx/gui/plot3d/scene/utils.py')
| -rw-r--r-- | src/silx/gui/plot3d/scene/utils.py | 662 |
1 files changed, 662 insertions, 0 deletions
diff --git a/src/silx/gui/plot3d/scene/utils.py b/src/silx/gui/plot3d/scene/utils.py new file mode 100644 index 0000000..c6cd129 --- /dev/null +++ b/src/silx/gui/plot3d/scene/utils.py @@ -0,0 +1,662 @@ +# 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 ################################################################# + +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 + + +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 + + +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) + + +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 + + +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 ######################################################################## + +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 + + +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 + + +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 ################################################################ + +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 + + +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 [] + + +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 + + +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 + + +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 ####################################################################### + +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() + + 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.) + + def move(self, step): + """Move the plane of step along the normal.""" + self.point += step * self.normal + + 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) |