Source code for fsleyes.gl.routines

#
# routines.py - A collection of disparate utility functions related to OpenGL.
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module contains a collection of miscellaneous OpenGL and geometric
routines.
"""


from __future__ import division

import                    logging
import                    contextlib
import collections.abc as abc
import itertools       as it

import OpenGL.GL       as gl
import numpy           as np

import fsl.transform.affine as affine


log = logging.getLogger(__name__)


[docs]def clear(bgColour): """Clears the current frame buffer, and does some standard setup operations. """ # set the background colour gl.glClearColor(*bgColour) # clear the buffer gl.glClear(gl.GL_COLOR_BUFFER_BIT | gl.GL_DEPTH_BUFFER_BIT) # enable transparency gl.glEnable(gl.GL_BLEND) gl.glBlendFunc(gl.GL_SRC_ALPHA, gl.GL_ONE_MINUS_SRC_ALPHA)
[docs]@contextlib.contextmanager def enabled(capabilities, enable=True): """This function can be used as a context manager to temporarily enable/disable one or more GL capabilities (i.e. something that you pass to ``glEnable`` and ``glDisable``, or ``glEnableClientState`` and ``glDisableClientState``), for a piece of code, and restore their previous state afterwards. :arg capabilities: One or more OpenGL capabilities to be temporarily enabled/disabled, e.g. ``GL_BLEND``, ``GL_TEXTURE_2D``, etc. :arg enable: Whether the capabilities are to be enabled (the default) or disabled. """ # Capabilities which are used via # glEnableClientState/glDisableClientState, # rather than glEnable/glDisable clientCapabilities = [ gl.GL_VERTEX_ARRAY, gl.GL_NORMAL_ARRAY, gl.GL_COLOR_ARRAY, gl.GL_SECONDARY_COLOR_ARRAY, gl.GL_EDGE_FLAG_ARRAY, gl.GL_INDEX_ARRAY, gl.GL_FOG_COORD_ARRAY, gl.GL_TEXTURE_COORD_ARRAY] if not isinstance(capabilities, abc.Sequence): capabilities = [capabilities] # Build lists of pre/post-yield # functions, one pre and post # for each capability, being # one of: # # - glEnable # - glDisable # - glEnableClientState # - glDisableClientState # - noop (if the capability's # current state is already # in the requested state) def noop(*a): pass pres = [] posts = [] for c in capabilities: if bool(gl.glIsEnabled(c)) == enable: pre = noop post = noop else: if c in clientCapabilities: pre = gl.glEnableClientState post = gl.glDisableClientState else: pre = gl.glEnable post = gl.glDisable if not enable: pre, post = post, pre pres .append(pre) posts.append(post) [p(c) for p, c in zip(pres, capabilities)] yield [p(c) for p, c in zip(posts, capabilities)]
[docs]@contextlib.contextmanager def disabled(capabilities): """Can be used as a context manager to temporarily disable the given GL capabilities - see :func:`enabled`. """ with enabled(capabilities, False): yield
[docs]def show2D(xax, yax, width, height, lo, hi, flipx=False, flipy=False): """Configures the OpenGL viewport for 2D othorgraphic display. :arg xax: Index (into ``lo`` and ``hi``) of the axis which corresponds to the horizontal screen axis. :arg yax: Index (into ``lo`` and ``hi``) of the axis which corresponds to the vertical screen axis. :arg width: Canvas width in pixels. :arg height: Canvas height in pixels. :arg lo: Tuple containing the mininum ``(x, y, z)`` display coordinates. :arg hi: Tuple containing the maxinum ``(x, y, z)`` display coordinates. :arg flipx: If ``True``, the x axis is inverted. :arg flipy: If ``True``, the y axis is inverted. """ zax = 3 - xax - yax xmin, xmax = lo[xax], hi[xax] ymin, ymax = lo[yax], hi[yax] zmin, zmax = lo[zax], hi[zax] projmat = np.eye(4, dtype=np.float32) if flipx: projmat[0, 0] = -1 if flipy: projmat[1, 1] = -1 gl.glViewport(0, 0, width, height) gl.glMatrixMode(gl.GL_PROJECTION) gl.glLoadMatrixf(projmat) zdist = max(abs(zmin), abs(zmax)) log.debug('Configuring orthographic viewport: ' 'X: [{} - {}] Y: [{} - {}] Z: [{} - {}]'.format( xmin, xmax, ymin, ymax, -zdist, zdist)) gl.glOrtho(xmin, xmax, ymin, ymax, -zdist, zdist) gl.glMatrixMode(gl.GL_MODELVIEW) gl.glLoadIdentity() # Rotate world space so the displayed slice # is visible and correctly oriented # TODO There's got to be a more generic way # to perform this rotation. This will break # if I add functionality allowing the user # to specifty the x/y axes on initialisation. if zax == 0: gl.glRotatef(270, 1, 0, 0) gl.glRotatef(270, 0, 0, 1) elif zax == 1: gl.glRotatef(270, 1, 0, 0)
[docs]def lookAt(eye, centre, up): """Replacement for ``gluLookAt`. Creates a transformation matrix which transforms the display coordinate system such that a camera at position (0, 0, 0), and looking towards (0, 0, -1), will see a scene as if from position ``eye``, oriented ``up``, and looking towards ``centre``. See: https://www.khronos.org/registry/OpenGL-Refpages/gl2.1/xhtml/gluLookAt.xml """ eye = np.array(eye) centre = np.array(centre) up = np.array(up) proj = np.eye(4) forward = centre - eye forward /= np.sqrt(np.dot(forward, forward)) right = np.cross(forward, up) right /= np.sqrt(np.dot(right, right)) up = np.cross(right, forward) up /= np.sqrt(np.dot(up, up)) proj[0, :3] = right proj[1, :3] = up proj[2, :3] = -forward eye = affine.scaleOffsetXform(1, -eye) proj = affine.concat(proj, eye) return proj
[docs]def ortho(lo, hi, width, height, zoom): """Generates an orthographic projection matrix. The display coordinate system origin ``(0, 0, 0)`` is mapped to the centre of the clipping space. - The horizontal axis is scaled to:: [-(hi[0] - lo[0]) / 2, (hi[0] - lo[0]) / 2] - The vertical axis is scaled to:: [-(hi[1] - lo[1]) / 2, (hi[1] - lo[1]) / 2] :arg lo: Low ``(x, y, z)`` bounds. :arg hi: High ``(x, y, z)`` bounds. :arg width: Canvas width in pixels :arg height: Canvas height in pixels :arg zoom: Zoom factor. Required to determine suitable near and far clipping plane locations. :returns: A tuple containing: - The ``(4, 4)`` projection matrix - A list of three ``(min, max)`` tuples, defining the limits for each axis. """ lo = np.array(lo, copy=False) hi = np.array(hi, copy=False) xmin, ymin = lo[:2] xmax, ymax = hi[:2] zmin, zmax = min(lo), max(hi) xmin, xmax, ymin, ymax = preserveAspectRatio( width, height, xmin, xmax, ymin, ymax) xlen = xmax - xmin ylen = ymax - ymin zlen = np.sqrt(np.sum((hi - lo) ** 2)) * zoom xhalf = xlen / 2.0 yhalf = ylen / 2.0 xmax = xhalf xmin = -xhalf ymax = yhalf ymin = -yhalf zmin = -zlen zmax = +zlen projmat = np.eye(4, dtype=np.float32) projmat[0, 0] = 2 / (xmax - xmin) projmat[1, 1] = 2 / (ymax - ymin) projmat[2, 2] = -2 / (zmax - zmin) projmat[0, 3] = -(xmax + xmin) / (xmax - xmin) projmat[1, 3] = -(ymax + ymin) / (ymax - ymin) projmat[2, 3] = -(zmax + zmin) / (zmax - zmin) limits = [(xmin, xmax), (ymin, ymax), (zmin, zmax)] return projmat, limits
[docs]def adjust(x, y, w, h): """Adjust the given ``x`` and ``y`` values by the aspect ratio defined by the given ``w`` and ``h`` values. """ xmin, xmax, ymin, ymax = preserveAspectRatio(w, h, 0, x, 0, y) return (xmax - xmin), (ymax - ymin)
[docs]def preserveAspectRatio(width, height, xmin, xmax, ymin, ymax, grow=True): """Adjusts the given x/y limits so that they can be displayed on a display of the given ``width`` and ``height``, while preserving the aspect ratio. :arg width: Display width :arg height: Display height :arg xmin: Low x limit :arg xmax: High x limit :arg ymin: Low y limit :arg ymax: High y limit :arg grow: If ``True`` (the default), the x/y limits are expanded to preserve the aspect ratio. Otherwise, they are shrunken. """ xlen = xmax - xmin ylen = ymax - ymin if np.any(np.isclose((width, height, xlen, ylen), 0)): return xmin, xmax, ymin, ymax # These ratios are used to determine whether # we need to expand the display range to # preserve the image aspect ratio. dispRatio = xlen / ylen canvasRatio = float(width) / height # the canvas is too wide - we need # to grow/shrink the display width, # thus effectively shrinking/growing # the display along the horizontal # axis if canvasRatio > dispRatio: newxlen = width * (ylen / height) if grow: offset = 0.5 * (newxlen - xlen) else: offset = -0.5 * (newxlen - xlen) xmin = xmin - offset xmax = xmax + offset # the canvas is too high - we need # to expand/shrink the display height elif canvasRatio < dispRatio: newylen = height * (xlen / width) if grow: offset = 0.5 * (newylen - ylen) else: offset = -0.5 * (newylen - ylen) ymin = ymin - offset ymax = ymax + offset return xmin, xmax, ymin, ymax
[docs]def pointGrid(shape, resolution, xform, xax, yax, origin='centre', bbox=None): """Calculates a uniform grid of points, in the display coordinate system (as specified by the given :class:`.Display` object properties) along the x-y plane (as specified by the xax/yax indices), at which the given image should be sampled for display purposes. This function returns a tuple containing: - a numpy array of shape ``(N, 3)``, containing the coordinates of the centre of every sampling point in the display coordinate system. - the horizontal distance (along xax) between adjacent points - the vertical distance (along yax) between adjacent points - The number of samples along the horizontal axis (xax) - The number of samples along the vertical axis (yax) :arg shape: The shape of the data to be sampled. :arg resolution: The desired resolution in display coordinates, along each display axis. :arg xform: A transformation matrix which converts from data coordinates to the display coordinate system. :arg xax: The horizontal display coordinate system axis (0, 1, or 2). :arg yax: The vertical display coordinate system axis (0, 1, or 2). :arg origin: ``centre`` or ``corner``. See the :func:`.affine.axisBounds` function. :arg bbox: An optional sequence of three ``(low, high)`` values, defining the bounding box in the display coordinate system which should be considered - the generated grid will be constrained to lie within this bounding box. """ xres = resolution[xax] yres = resolution[yax] # These values give the min/max x/y # values of a bounding box which # encapsulates the entire image, # in the display coordinate system xmin, xmax = affine.axisBounds(shape, xform, xax, origin, boundary=None) ymin, ymax = affine.axisBounds(shape, xform, yax, origin, boundary=None) # Number of samples along each display # axis, given the requested resolution xNumSamples = int(np.floor((xmax - xmin) / xres)) yNumSamples = int(np.floor((ymax - ymin) / yres)) # adjust the x/y resolution so # the samples fit exactly into # the data bounding box xres = (xmax - xmin) / xNumSamples yres = (ymax - ymin) / yNumSamples # Calculate the locations of every # sample point in display space worldX = np.linspace(xmin + 0.5 * xres, xmax - 0.5 * xres, xNumSamples) worldY = np.linspace(ymin + 0.5 * yres, ymax - 0.5 * yres, yNumSamples) # Apply bounding box constraint # if it has been provided if bbox is not None: xoff = 0.5 * xres yoff = 0.5 * yres xmin = max((xmin, bbox[xax][0] - xoff)) xmax = min((xmax, bbox[xax][1] + xoff)) ymin = max((ymin, bbox[yax][0] - yoff)) ymax = min((ymax, bbox[yax][1] + yoff)) worldX = worldX[(worldX >= xmin) & (worldX <= xmax)] worldY = worldY[(worldY >= ymin) & (worldY <= ymax)] # Generate the coordinates worldX, worldY = np.meshgrid(worldX, worldY) # reshape them to N*3 coords = np.zeros((worldX.size, 3), dtype=np.float32) coords[:, xax] = worldX.flatten() coords[:, yax] = worldY.flatten() return coords, xres, yres, xNumSamples, yNumSamples
[docs]def pointGrid3D(shape, xform=None, origin='centre', bbox=None): """Generates a 3D grid of points into an image of the given ``shape``, with the given ``xform`` defining the index to display coordinate transform. note: Not implemented properly yet. """ coords = np.indices(shape).transpose(1, 2, 3, 0).reshape(-1, 3) if xform is not None: coords = affine.transform(coords, xform) return coords
[docs]def samplePointsToTriangleStrip(coords, xpixdim, ypixdim, xlen, ylen, xax, yax): """Given a regular 2D grid of points at which an image is to be sampled (for example, that generated by the :func:`pointGrid` function above), converts those points into an OpenGL vertex triangle strip. A grid of ``M*N`` points is represented by ``M*2*(N + 1)`` vertices. For example, this image represents a 4*3 grid, with periods representing vertex locations:: .___.___.___.___. | | | | | | | | | | .---.---.---.---. .___.___.__ .___. | | | | | | | | | | .---.---.---.---. .___.___.___.___. | | | | | | | | | | .___.___.___.___. Vertex locations which are vertically adjacent represent the same point in space. Such vertex pairs are unable to be combined because, in OpenGL, they must be represented by distinct vertices (we can't apply multiple colours/texture coordinates to a single vertex location) So we have to repeat these vertices in order to achieve accurate colouring of each voxel. We draw each horizontal row of samples one by one, using two triangles to draw each voxel. In order to eliminate the need to specify six vertices for every voxel, and hence to reduce the amount of memory used, we are using a triangle strip to draw each row of voxels. This image depicts a triangle strip used to draw a row of three samples (periods represent vertex locations):: 1 3 5 7 . . . . |\\ |\\ |\\ | | \\| \\| \\| . . . . 0 2 4 6 In order to use a single OpenGL call to draw multiple non-contiguous voxel rows, between every column we add a couple of 'dummy' vertices, which will then be interpreted by OpenGL as 'degenerate triangles', and will not be drawn. So in reality, a 4*3 slice would be drawn as follows (with vertices labelled from ``[a-z0-9]``:: v x z 1 33 |\\ |\\ |\\ |\\ | | \\| \\| \\| \\| uu w y 0 2 l n p r tt |\\ |\\ |\\ |\\ | | \\| \\| \\| \\| kk m o q s b d f h jj |\\ |\\ |\\ |\\ | | \\| \\| \\| \\| a c e g i These repeated/degenerate vertices are dealt with by using a vertex index array. See these links for good overviews of triangle strips and degenerate triangles in OpenGL: - http://www.learnopengles.com/tag/degenerate-triangles/ - http://en.wikipedia.org/wiki/Triangle_strip A tuple is returned containing: - A 2D ``numpy.float32`` array of shape ``(2 * (xlen + 1) * ylen), 3)``, containing the coordinates of all triangle strip vertices which represent the entire grid of sample points. - A 2D ``numpy.float32`` array of shape ``(2 * (xlen + 1) * ylen), 3)``, containing the centre of every grid, to be used for texture coordinates/colour lookup. - A 1D ``numpy.uint32`` array of size ``ylen * (2 * (xlen + 1) + 2) - 2`` containing indices into the first array, defining the order in which the vertices need to be rendered. There are more indices than vertex coordinates due to the inclusion of repeated/degenerate vertices. :arg coords: N*3 array of points, the sampling locations. :arg xpixdim: Length of one sample along the horizontal axis. :arg ypixdim: Length of one sample along the vertical axis. :arg xlen: Number of samples along the horizontal axis. :arg ylen: Number of samples along the vertical axis. :arg xax: Display coordinate system axis which corresponds to the horizontal screen axis. :arg yax: Display coordinate system axis which corresponds to the vertical screen axis. """ coords = coords.reshape(ylen, xlen, 3) xlen = int(xlen) ylen = int(ylen) # Duplicate every row - each voxel # is defined by two vertices coords = coords.repeat(2, 0) texCoords = np.array(coords, dtype=np.float32) worldCoords = np.array(coords, dtype=np.float32) # Add an extra column at the end # of the world coordinates worldCoords = np.append(worldCoords, worldCoords[:, -1:, :], 1) worldCoords[:, -1, xax] += xpixdim # Add an extra column at the start # of the texture coordinates texCoords = np.append(texCoords[:, :1, :], texCoords, 1) # Move the x/y world coordinates to the # sampling point corners (the texture # coordinates remain in the voxel centres) worldCoords[ :, :, xax] -= 0.5 * xpixdim worldCoords[ ::2, :, yax] -= 0.5 * ypixdim worldCoords[1::2, :, yax] += 0.5 * ypixdim vertsPerRow = 2 * (xlen + 1) dVertsPerRow = 2 * (xlen + 1) + 2 nindices = ylen * dVertsPerRow - 2 indices = np.zeros(nindices, dtype=np.uint32) for yi, xi in it.product(range(ylen), range(xlen + 1)): ii = yi * dVertsPerRow + 2 * xi vi = yi * vertsPerRow + xi indices[ii] = vi indices[ii + 1] = vi + xlen + 1 # add degenerate vertices at the end # every row (but not needed for last # row) if xi == xlen and yi < ylen - 1: indices[ii + 2] = vi + xlen + 1 indices[ii + 3] = (yi + 1) * vertsPerRow worldCoords = worldCoords.reshape((xlen + 1) * (2 * ylen), 3) texCoords = texCoords .reshape((xlen + 1) * (2 * ylen), 3) return worldCoords, texCoords, indices
[docs]def voxelGrid(points, xax, yax, xpixdim, ypixdim): """Given a ``N*3`` array of ``points`` (assumed to be voxel coordinates), creates an array of vertices which can be used to render each point as an unfilled rectangle. :arg points: An ``N*3`` array of voxel xyz coordinates :arg xax: XYZ axis index that maps to the horizontal scren axis :arg yax: XYZ axis index that maps to the vertical scren axis :arg xpixdim: Length of a voxel along the x axis. :arg ypixdim: Length of a voxel along the y axis. """ if len(points.shape) == 1: points = points.reshape(1, 3) npoints = points.shape[0] vertices = np.repeat(np.array(points, dtype=np.float32), 4, axis=0) xpixdim = xpixdim / 2.0 ypixdim = ypixdim / 2.0 # bottom left corner vertices[ ::4, xax] -= xpixdim vertices[ ::4, yax] -= ypixdim # bottom right vertices[1::4, xax] += xpixdim vertices[1::4, yax] -= ypixdim # top left vertices[2::4, xax] -= xpixdim vertices[2::4, yax] += ypixdim # top right vertices[3::4, xax] += xpixdim vertices[3::4, yax] += ypixdim # each square is rendered as four lines indices = np.array([0, 1, 0, 2, 1, 3, 2, 3], dtype=np.uint32) indices = np.tile(indices, npoints) indices = (indices.T + np.repeat(np.arange(0, npoints * 4, 4, dtype=np.uint32), 8)).T return vertices, indices
[docs]def voxelBlock(*args, **kwargs): """Generates a ``numpy`` array containing all ones, centered at the specified voxel. :arg dtype: The data type of the returned ``numpy`` array. Defaults to ``uint8``. All other arguments are passed through to the :func:`voxelBox` function - see that function for details on the arguments. :returns: A tuple containing: - The ``numpy`` array - Voxel coordinates specifying the offset of the position of this array in the image. """ dtype = kwargs.pop('dtype', np.uint8) corners = voxelBox(*args, **kwargs) los = corners.min(axis=0) his = corners.max(axis=0) lens = his - los block = np.ones(lens, dtype=dtype) return block, los
[docs]def voxelBox(voxel, shape, dims, boxSize, axes=(0, 1, 2), bias=None, bounded=True): """Generates a 'box', a cuboid region, in a voxel coordinate system, centered at a specific voxel. The corners of the box are returnd as a ``numpy`` array of shape ``(8, 3)``. :arg voxel: Coordinates of the voxel around which the block is to be centred. :arg shape: Shape of the image in which the block is to be located. :arg dims: Size of the image voxels along each dimension. :arg boxSize: Desired width/height/depth of the box in scaled voxels. May be either a single value, or a sequence of three values. The box size is forced to be at least one voxel in each dimension. :arg axes: Axes along which the block is to be located. :arg bias: ``low``, ``high``, or ``None``. The box can only be centered on the ``voxel`` if the specified ``boxSize`` results in an odd number of voxels along each voxel axis. If this is not the case, more voxels must be added to one side of the box centre. You can specify that these voxels are added to the ``high`` side (i.e. larger voxel coordinate) or the ``low`` side of the voxel. Specifying ``None`` will force the box to have an odd number of voxels along each axis, and thus have the ``voxel`` in the true centre of the box. :arg bounded: If ``True`` (the default), and the specified voxel would result in part of the block being located outside of the image shape, the block is truncated to fit inside the image bounds. :returns: A ``(8, 3)`` numpy array containing the corners of the box, or ``None`` if a box of the requested size/shape cannot be defined in the specifiied coordinate system. """ if not isinstance(boxSize, abc.Iterable): boxSize = [boxSize] * 3 # Limit the box span to the specified # axes (e.g. a 2D plane) for i in range(3): if i not in axes: boxSize[i] = dims[i] voxel = np.array(voxel) dims = np.array(dims[:3]) shape = np.array(shape[:3]) boxSize = np.array(boxSize[:3]) # The block size must be at least # one voxel across each dimension boxSize = np.clip(boxSize, dims, None) # Voxel location and box low/ # high bounds in scaled voxels. # # Note that we are assuming that # voxel coordinates correspond to # the voxel centre here # (voxel + 0.5). This has the effect # that the returned box vertices will # be integer voxel coordinates (i.e. # on voxel borders). scaledVoxel = (voxel + 0.5) * dims scaledLo = scaledVoxel - boxSize / 2.0 scaledHi = scaledVoxel + boxSize / 2.0 # Scale the low/high bounds back # into voxel coordinates, and # round them up or down according # to the bias setting. if bias == 'low': lo = np.floor(scaledLo / dims) hi = np.floor(scaledHi / dims) elif bias == 'high': lo = np.ceil(scaledLo / dims) hi = np.ceil(scaledHi / dims) else: lo = np.floor(scaledLo / dims) hi = np.ceil( scaledHi / dims) # Crop the box to the # image space if necessary if bounded: lo = np.maximum(lo, 0) hi = np.minimum(hi, shape) # requesting a box that is smaller than # one voxel along at least one axis if np.any(hi <= lo): return None box = np.zeros((8, 3), dtype=np.uint32) box[0, :] = lo[0], lo[1], lo[2] box[1, :] = lo[0], lo[1], hi[2] box[2, :] = lo[0], hi[1], lo[2] box[3, :] = lo[0], hi[1], hi[2] box[4, :] = hi[0], lo[1], lo[2] box[5, :] = hi[0], lo[1], hi[2] box[6, :] = hi[0], hi[1], lo[2] box[7, :] = hi[0], hi[1], hi[2] return box
[docs]def slice2D(dataShape, xax, yax, zpos, voxToDisplayMat, displayToVoxMat, geometry='triangles', origin='centre', bbox=None): """Generates and returns vertices which denote a slice through an array of the given ``dataShape``, parallel to the plane defined by the given ``xax`` and ``yax`` and at the given z position, in the space defined by the given ``voxToDisplayMat``. If ``geometry`` is ``triangles`` (the default), six vertices are returned, arranged as follows:: 4----5 1 \\ | |\\ \\ | | \\ \\| | \\ 3 0----2 Otherwise, if geometry is ``square``, four vertices are returned, arranged as follows:: 3---2 | | | | | | 0---1 If ``origin`` is set to ``centre`` (the default), it is assumed that a voxel at location ``(x, y, z)`` is located in the space:: (x - 0.5 : x + 0.5, y - 0.5 : y + 0.5, z - 0.5 : z + 0.5) Otherwise, if ``origin`` is set to ``corner``, a voxel at location ``(x, y, z)`` is assumed to be located in the space:: (x : x + 1, y : y + 1, z : z + 1) :arg dataShape: Number of elements along each dimension in the image data. :arg xax: Index of display axis which corresponds to the horizontal screen axis. :arg yax: Index of display axis which corresponds to the vertical screen axis. :arg zpos: Position of the slice along the screen z axis. :arg voxToDisplayMat: Affine transformation matrix which transforms from voxel/array indices into the display coordinate system. :arg displayToVoxMat: Inverse of the ``voxToDisplayMat``. :arg geometry: ``square`` or ``triangle``. :arg origin: ``centre`` or ``corner``. See the :func:`.affine.axisBounds` function. :arg bbox: An optional sequence of three ``(low, high)`` values, defining the bounding box in the display coordinate system which should be considered - the generated grid will be constrained to lie within this bounding box. Returns a tuple containing: - A ``N*3`` ``numpy.float32`` array containing the vertex locations of a slice through the data, where ``N=6`` if ``geometry=triangles``, or ``N=4`` if ``geometry=square``, - A ``N*3`` ``numpy.float32`` array containing the voxel coordinates that correspond to the vertex locations. """ zax = 3 - xax - yax xmin, xmax = affine.axisBounds( dataShape, voxToDisplayMat, xax, origin, boundary=None) ymin, ymax = affine.axisBounds( dataShape, voxToDisplayMat, yax, origin, boundary=None) if bbox is not None: bbxmin = bbox[xax][0] bbxmax = bbox[xax][1] bbymin = bbox[yax][0] bbymax = bbox[yax][1] # The returned vertices and voxCoords # have to be aligned, so we need to # clamp the bounding box limits to the # nearest voxel boundary to preserve # this alignment. # Voxel lengths along x/y axes xvlen = (xmax - xmin) / dataShape[xax] yvlen = (ymax - ymin) / dataShape[yax] # Clamp the bbox limits to the # nearest voxel boundaries bbxmin = xmin + np.floor((bbxmin - xmin) / xvlen) * xvlen bbxmax = xmin + np.ceil( (bbxmax - xmin) / xvlen) * xvlen bbymin = ymin + np.floor((bbymin - ymin) / yvlen) * yvlen bbymax = ymin + np.ceil( (bbymax - ymin) / yvlen) * yvlen xmin = max((xmin, bbxmin)) xmax = min((xmax, bbxmax)) ymin = max((ymin, bbymin)) ymax = min((ymax, bbymax)) if geometry == 'triangles': vertices = np.zeros((6, 3), dtype=np.float32) vertices[ 0, [xax, yax]] = [xmin, ymin] vertices[ 1, [xax, yax]] = [xmax, ymin] vertices[ 2, [xax, yax]] = [xmin, ymax] vertices[ 3, [xax, yax]] = [xmax, ymin] vertices[ 4, [xax, yax]] = [xmax, ymax] vertices[ 5, [xax, yax]] = [xmin, ymax] elif geometry == 'square': vertices = np.zeros((4, 3), dtype=np.float32) vertices[ 0, [xax, yax]] = [xmin, ymin] vertices[ 1, [xax, yax]] = [xmax, ymin] vertices[ 2, [xax, yax]] = [xmax, ymax] vertices[ 3, [xax, yax]] = [xmin, ymax] else: raise ValueError('Unrecognised geometry type: {}'.format(geometry)) vertices[:, zax] = zpos voxCoords = affine.transform(vertices, displayToVoxMat) return vertices, voxCoords
[docs]def boundingBox(dataShape, voxToDisplayMat, displayToVoxMat, geometry='triangles', origin='centre', bbox=None): """Generates a bounding box to represent a 3D image of the given shape, in the coordinate system defined by the ``voxToDisplayMat`` affine. See the :func:`slice2D` function for details on the arguments. Returns a tuple containing: - A ``N*3`` ``numpy.float32`` array containing the vertex locations of a bounding box ``N=36`` if ``geometry=triangles``, or ``N=24`` if ``geometry=square``, - A ``N*3`` ``numpy.float32`` array containing the voxel coordinates that correspond to the vertex locations. """ xlo, ylo, zlo = (0, 0, 0) xhi, yhi, zhi = dataShape if origin == 'centre': xlo, ylo, zlo = xlo - 0.5, ylo - 0.5, zlo - 0.5 xhi, yhi, zhi = xhi - 0.5, yhi - 0.5, zhi - 0.5 # If the voxel -> display transformation # matrix contains negative scales, we need # to invert the voxel coordinate ranges to # ensure a correct triangle winding order # (so that back faces can be culled). scales = affine.decompose(voxToDisplayMat)[0] if scales[0] < 0: xlo, xhi = xhi, xlo if scales[1] < 0: ylo, yhi = yhi, ylo if scales[2] < 0: zlo, zhi = zhi, zlo if geometry == 'triangles': voxCoords = np.zeros((36, 3), dtype=np.float32) voxCoords[ 0, :] = (xlo, yhi, zlo) voxCoords[ 1, :] = (xhi, ylo, zlo) voxCoords[ 2, :] = (xlo, ylo, zlo) voxCoords[ 3, :] = (xlo, yhi, zlo) voxCoords[ 4, :] = (xhi, yhi, zlo) voxCoords[ 5, :] = (xhi, ylo, zlo) voxCoords[ 6, :] = (xlo, ylo, zhi) voxCoords[ 7, :] = (xhi, ylo, zhi) voxCoords[ 8, :] = (xlo, yhi, zhi) voxCoords[ 9, :] = (xlo, yhi, zhi) voxCoords[10, :] = (xhi, ylo, zhi) voxCoords[11, :] = (xhi, yhi, zhi) voxCoords[12, :] = (xlo, ylo, zlo) voxCoords[13, :] = (xhi, ylo, zlo) voxCoords[14, :] = (xlo, ylo, zhi) voxCoords[15, :] = (xlo, ylo, zhi) voxCoords[16, :] = (xhi, ylo, zlo) voxCoords[17, :] = (xhi, ylo, zhi) voxCoords[18, :] = (xlo, yhi, zlo) voxCoords[19, :] = (xhi, yhi, zlo) voxCoords[20, :] = (xlo, yhi, zhi) voxCoords[21, :] = (xlo, yhi, zhi) voxCoords[22, :] = (xhi, yhi, zlo) voxCoords[23, :] = (xhi, yhi, zhi) voxCoords[18, :] = (xlo, yhi, zhi) voxCoords[19, :] = (xhi, yhi, zlo) voxCoords[20, :] = (xlo, yhi, zlo) voxCoords[21, :] = (xlo, yhi, zhi) voxCoords[22, :] = (xhi, yhi, zhi) voxCoords[23, :] = (xhi, yhi, zlo) voxCoords[24, :] = (xlo, ylo, zhi) voxCoords[25, :] = (xlo, yhi, zlo) voxCoords[26, :] = (xlo, ylo, zlo) voxCoords[27, :] = (xlo, ylo, zhi) voxCoords[28, :] = (xlo, yhi, zhi) voxCoords[29, :] = (xlo, yhi, zlo) voxCoords[30, :] = (xhi, ylo, zlo) voxCoords[31, :] = (xhi, yhi, zlo) voxCoords[32, :] = (xhi, ylo, zhi) voxCoords[33, :] = (xhi, ylo, zhi) voxCoords[34, :] = (xhi, yhi, zlo) voxCoords[35, :] = (xhi, yhi, zhi) elif geometry == 'square': voxCoords = np.zeros((24, 3), dtype=np.float32) voxCoords[ 0, :] = (xlo, ylo, zlo) voxCoords[ 1, :] = (xhi, ylo, zlo) voxCoords[ 2, :] = (xhi, yhi, zlo) voxCoords[ 3, :] = (xlo, yhi, zlo) voxCoords[ 4, :] = (xlo, ylo, zhi) voxCoords[ 5, :] = (xhi, ylo, zhi) voxCoords[ 6, :] = (xhi, yhi, zhi) voxCoords[ 7, :] = (xlo, yhi, zhi) voxCoords[ 8, :] = (xlo, ylo, zlo) voxCoords[ 9, :] = (xhi, ylo, zlo) voxCoords[10, :] = (xhi, ylo, zhi) voxCoords[11, :] = (xlo, ylo, zhi) voxCoords[12, :] = (xlo, yhi, zlo) voxCoords[13, :] = (xhi, yhi, zlo) voxCoords[14, :] = (xhi, yhi, zhi) voxCoords[15, :] = (xlo, yhi, zhi) voxCoords[16, :] = (xlo, ylo, zlo) voxCoords[17, :] = (xlo, yhi, zlo) voxCoords[18, :] = (xlo, yhi, zhi) voxCoords[19, :] = (xlo, ylo, zhi) voxCoords[20, :] = (xhi, ylo, zlo) voxCoords[21, :] = (xhi, yhi, zlo) voxCoords[22, :] = (xhi, yhi, zhi) voxCoords[23, :] = (xhi, ylo, zhi) else: raise ValueError('Unrecognised geometry type: {}'.format(geometry)) vertices = affine.transform(voxCoords, voxToDisplayMat) return vertices, voxCoords
[docs]def subsample(data, resolution, pixdim=None, volume=None): """Samples the given 3D data according to the given resolution. Returns a tuple containing: - A 3D numpy array containing the sub-sampled data. - A tuple containing the ``(x, y, z)`` starting indices of the sampled data. - A tuple containing the ``(x, y, z)`` steps of the sampled data. :arg data: The data to be sampled. :arg resolution: Sampling resolution, proportional to the values in ``pixdim``. :arg pixdim: Length of each dimension in the input data (defaults to ``(1.0, 1.0, 1.0)``). :arg volume: If the image is a 4D volume, the volume index of the 3D image to be sampled. """ if pixdim is None: pixdim = (1.0, 1.0, 1.0) if volume is None: volume = slice(None, None, None) xstep = int(np.round(resolution / pixdim[0])) ystep = int(np.round(resolution / pixdim[1])) zstep = int(np.round(resolution / pixdim[2])) if xstep < 1: xstep = 1 if ystep < 1: ystep = 1 if zstep < 1: zstep = 1 xstart = int(np.floor((xstep - 1) / 2)) ystart = int(np.floor((ystep - 1) / 2)) zstart = int(np.floor((zstep - 1) / 2)) if xstart >= data.shape[0]: xstart = data.shape[0] - 1 if ystart >= data.shape[1]: ystart = data.shape[1] - 1 if zstart >= data.shape[2]: zstart = data.shape[2] - 1 if len(data.shape) > 3: sample = data[xstart::xstep, ystart::ystep, zstart::zstep, volume] else: sample = data[xstart::xstep, ystart::ystep, zstart::zstep] return sample, (xstart, ystart, zstart), (xstep, ystep, zstep)
[docs]def broadcast(vertices, indices, zposes, xforms, zax): """Given a set of vertices and indices (assumed to be 2D representations of some geometry in a 3D space, with the depth axis specified by ``zax``), replicates them across all of the specified Z positions, applying the corresponding transformation to each set of vertices. :arg vertices: Vertex array (a ``N*3`` numpy array). :arg indices: Index array. :arg zposes: Positions along the depth axis at which the vertices are to be replicated. :arg xforms: Sequence of transformation matrices, one for each Z position. :arg zax: Index of the 'depth' axis Returns three values: - A numpy array containing all of the generated vertices - A numpy array containing the original vertices for each of the generated vertices, which may be used as texture coordinates - A new numpy array containing all of the generated indices. """ vertices = np.array(vertices) indices = np.array(indices) nverts = vertices.shape[0] nidxs = indices.shape[ 0] allTexCoords = np.zeros((nverts * len(zposes), 3), dtype=np.float32) allVertCoords = np.zeros((nverts * len(zposes), 3), dtype=np.float32) allIndices = np.zeros( nidxs * len(zposes), dtype=np.uint32) for i, (zpos, xform) in enumerate(zip(zposes, xforms)): vertices[:, zax] = zpos vStart = i * nverts vEnd = vStart + nverts iStart = i * nidxs iEnd = iStart + nidxs allIndices[ iStart:iEnd] = indices + i * nverts allTexCoords[ vStart:vEnd, :] = vertices allVertCoords[vStart:vEnd, :] = affine.transform(vertices, xform) return allVertCoords, allTexCoords, allIndices
[docs]def planeEquation(xyz1, xyz2, xyz3): """Calculates the equation of a plane which contains each of the given points. Returns a ``numpy`` array containing four values, the coefficients of the equation: :math:`a\\times x + b\\times y + c \\times z = d` for any point ``(x, y, z)`` that lies on the plane. See http://paulbourke.net/geometry/pointlineplane/ for details on plane equations. """ x1, y1, z1 = xyz1 x2, y2, z2 = xyz2 x3, y3, z3 = xyz3 eq = np.zeros(4, dtype=np.float64) eq[0] = (y1 * (z2 - z3)) + (y2 * (z3 - z1)) + (y3 * (z1 - z2)) eq[1] = (z1 * (x2 - x3)) + (z2 * (x3 - x1)) + (z3 * (x1 - x2)) eq[2] = (x1 * (y2 - y3)) + (x2 * (y3 - y1)) + (x3 * (y1 - y2)) eq[3] = -((x1 * ((y2 * z3) - (y3 * z2))) + (x2 * ((y3 * z1) - (y1 * z3))) + (x3 * ((y1 * z2) - (y2 * z1)))) return eq
[docs]def planeEquation2(origin, normal): """Calculates the equation of a plane equation from a normal vector and a single point on the plane. Returns a ``numpy`` array containing four values, the coefficients of the equation: See also :func:`planeEquation`. """ normal = affine.normalise(normal) ax, by, cz = np.array(origin) * normal eqn = np.zeros(4, dtype=np.float64) eqn[:3] = normal eqn[ 3] = -np.sum((ax, by, cz)) return eqn
[docs]def unitSphere(res): """Generates a unit sphere, as described in the *Sphere Generation* article, on Paul Bourke's excellent website: http://paulbourke.net/geometry/circlesphere/ :arg res: Resolution - the number of angles to sample. :returns: A tuple comprising: - a ``numpy.float32`` array of size ``(res**2, 3)`` containing a set of ``(x, y, z)`` vertices which define the ellipsoid surface. - A ``numpy.uint32`` array of size ``(4 * (res - 1)**2)`` containing a list of indices into the vertex array, defining a vertex ordering that can be used to draw the ellipsoid using the OpenGL ``GL_QUAD`` primitive type. .. todo:: Generate indices to use with ``GL_TRIANGLES`` instead of ``GL_QUADS``. """ # All angles to be sampled u = np.linspace(-np.pi / 2, np.pi / 2, res, dtype=np.float32) v = np.linspace(-np.pi, np.pi, res, dtype=np.float32) cosu = np.cos(u) cosv = np.cos(v) sinu = np.sin(u) sinv = np.sin(v) cucv = np.outer(cosu, cosv).T cusv = np.outer(cosu, sinv).T vertices = np.zeros((res ** 2, 3), dtype=np.float32) # All x coordinates are of the form cos(u) * cos(v), # y coordinates are of the form cos(u) * sin(v), # and z coordinates of the form sin(u). vertices[:, 0] = cucv.flatten() vertices[:, 1] = cusv.flatten() vertices[:, 2] = np.tile(sinu, res) # Generate a list of indices which join the # vertices so they can be used to draw the # sphere as GL_QUADs. # # The vertex locations for each quad follow # this pattern: # # 1. (u, v) # 2. (u + ustep, v) # 3. (u + ustep, v + vstep) # 4. (u, v + vstep) nquads = (res - 1) ** 2 quadIdxs = np.array([0, res, res + 1, 1], dtype=np.uint32) indices = np.tile(quadIdxs, nquads) indices += np.arange(nquads, dtype=np.uint32).repeat(4) indices += np.arange(res - 1, dtype=np.uint32).repeat(4 * (res - 1)) return vertices, indices
[docs]def fullUnitSphere(res): """Generates a unit sphere in the same way as :func:`unitSphere`, but returns all vertices, instead of the unique vertices and an index array. :arg res: Resolution - the number of angles to sample. :returns: A ``numpy.float32`` array of size ``(4 * (res - 1)**2, 3)`` containing the ``(x, y, z)`` vertices which can be used to draw a unit sphere (using the ``GL_QUADS`` primitive type). """ u = np.linspace(-np.pi / 2, np.pi / 2, res, dtype=np.float32) v = np.linspace(-np.pi, np.pi, res, dtype=np.float32) cosu = np.cos(u) cosv = np.cos(v) sinu = np.sin(u) sinv = np.sin(v) vertices = np.zeros(((res - 1) * (res - 1) * 4, 3), dtype=np.float32) cucv = np.outer(cosu[:-1], cosv[:-1]).flatten() cusv = np.outer(cosu[:-1], sinv[:-1]).flatten() cu1cv = np.outer(cosu[1:], cosv[:-1]).flatten() cu1sv = np.outer(cosu[1:], sinv[:-1]).flatten() cu1cv1 = np.outer(cosu[1:], cosv[1:]) .flatten() cu1sv1 = np.outer(cosu[1:], sinv[1:]) .flatten() cucv1 = np.outer(cosu[:-1], cosv[1:]) .flatten() cusv1 = np.outer(cosu[:-1], sinv[1:]) .flatten() su = np.repeat(sinu[:-1], res - 1) s1u = np.repeat(sinu[1:], res - 1) vertices.T[:, ::4] = [cucv, cusv, su] vertices.T[:, 1::4] = [cu1cv, cu1sv, s1u] vertices.T[:, 2::4] = [cu1cv1, cu1sv1, s1u] vertices.T[:, 3::4] = [cucv1, cusv1, su] return vertices
[docs]def unitCircle(res, triangles=False): """Generates ``res`` vertices which form a 2D circle, centered at (0, 0), and with radius 1. Returns the vertices as a ``numpy.float32`` array of shape ``(res, 2)`` If the ``triangles`` argument is ``True``, the vertices are generated with the assumption that they will be drawn as a ``GL_TRIANGLE_FAN``. """ step = (2 * np.pi) / res u = np.linspace(-np.pi, np.pi - step, res, dtype=np.float32) cosu = np.cos(u) sinu = np.sin(u) verts = np.vstack((sinu, cosu)).T if triangles: origin = np.zeros((1, 3), dtyp=np.float32) verts = np.concatenate((origin, verts)) return verts
[docs]def polygonIndices(nverts): """Generate triangle indices for simple 2D polygons. Given vertices describing a monotone polygon on a 2D plane, generates an index list into the vertices which can be used to draw the vertices as triangles. """ ntris = nverts - 2 nidxs = ntris * 3 indices = np.zeros(nidxs, dtype=np.uint32) indices[1] = 1 indices[2:-1] = np.repeat(np.arange(2, nverts - 1), 3) indices[-1] = nverts - 1 indices[::3] = 0 return indices