gyral_structure API

Describes white matter structure as a vector field. This vector field defines at every point the local fibre orientation and density.

• This vector field is defined as the sum of basis functions (defined in gyral_structure.basis).

• The coefficients of these basis functions are found by minimizing any cost function (see gyral_structure.cost).

• These cost functions generally require the vector field to be evaluated in certain voxels or vertices. These requests are handled by gyral_structure.request.

• Finally the algorithm used to compute this vector field is defined by gyral_structure.algorithm (e.g. cpu or gpu, pre-compute matrix or not).

• Once the best-fit coefficients have been found streamlines can be computed using gyral_structure.streamline.

You can view the code at https://git.fmrib.ox.ac.uk/ndcn0236/gyral_structure

• gyral_structure.basis: defining the basis functions
  • gyral_structure.basis.core module
  • gyral_structure.basis.evaluator module
  • gyral_structure.basis.grid module
  • gyral_structure.basis.radial module
  • gyral_structure.basis.tricubic module
  • gyral_structure.basis.trilinear module
  • gyral_structure.basis.unbounded module
• gyral_structure.cost package: defining the cost function
  • gyral_structure.cost.core module
  • gyral_structure.cost.dens module
  • gyral_structure.cost.meta module
  • gyral_structure.cost.orient module
  • gyral_structure.cost.param module
• Scripts included in gyral_structure
  • gs_mask: Creates a mask of the gyral white matter
  • gs_fit: Fits vector field inside gyral white matter
  • gs_deform_surface: Moves white/gray matter boundary to deep white matter
  • gs_dyad_vol: Evaluates vector field fibre density and orientation in a volume
  • gs_dyad_surf: Evaluates vector field fibre density and orientation on a surface
  • gs_track: Streamline tractography through vector field

gyral_structure.algorithm: how to compute the cost function

Use set() to alter how the cost functions is evaluated

gyral_structure.algorithm.full_chain(available=None)

Returns the ordered set of algorithms that should be tested for the basis function

Parameters

available – sequence of available basis functions in order of preference

Returns

re-ordered list of available basis functions based on preferences set in this module

class gyral_structure.algorithm.Algorithm

Bases: enum.Enum

Methods to compute the vector field based on the parameters.

matrix = 1

numpy = 2

numba = 3

matrix_cuda = 4

cuda = 5

convolve = 6

gyral_structure.algorithm.set(use_cuda: bool = None, store_matrix: bool = None, method: gyral_structure.algorithm.Algorithm = None, override=False)

Sets the default method to calculate the cost function.

Use in a with statement, for example

with set(use_cuda=False):
    pass

This causes any evaluators generated in the with-statement not to use the GPU.

Parameters

• use_cuda – if True prefers cuda-enabled methods, if False don’t use cuda (default: nothing is changed)

• store_matrix – if True, prefers to store the matrix prior to calculation, if False don’t precompute the matrix (default: nothing is changed)

• method – explicitly set the preferred algorithm (default: nothing is changed)

• override – if True overrides any previous settings (default: False)

class gyral_structure.algorithm.SettingsContext(use_cuda: bool = None, store_matrix: bool = None, method: gyral_structure.algorithm.Algorithm = None, override: bool = False)

Bases: object

__init__(use_cuda: bool = None, store_matrix: bool = None, method: gyral_structure.algorithm.Algorithm = None, override: bool = False)

Initialize self. See help(type(self)) for accurate signature.

gyral_structure.cuda: GPU support

Contains a large number of functions to evaluate the basis functions on the GPU

For pre-computed matrix multiplication a kernel to evaluate the vector field can be gained by providing the matrix to

• sparse_mult_prepare() for sparse matrices

• dense_mult_prepare() for dense matrices

On the fly computation of dense matrices can be implemented through CudaMatrixMult

gyral_structure.cuda.correct_type(arr)

If needed return converted array for transfer to GPU.

Arrays of the correct type will be returned unaltered. All float arrays will be made single/double precision (depending on cuda.dtype) All integer arrays will be made short (‘i4’)

Parameters

arr – input array

Returns

either arr itself if it is of the correct type or a converted arr to the correct type

gyral_structure.cuda.to_gpu_correct(arr)

Returns the GPUArray pointing to the input arr on the GPU

This algorithm ensures that the array has the correct type (converting if if necessary) and then moves the data to the GPu.

Parameters

arr – numpy data array

Returns

pointer to the data on the GPU

gyral_structure.cuda.sparse_mult_prepare(sparse)

Returns a function that does the GPU multiplication for the provided sparse matrix.

gyral_structure.cuda.dense_mult_prepare(dense)

Prepares a dense matrix multiplication.

Parameters

dense – (M, N) array

Returns

function that does the dense matrix multiplication

gyral_structure.cuda.wrap_get_element(code)

Defines kernel functions if a get_element CUDA function has been defined.

class gyral_structure.cuda.CudaMatrixMult(code, npos, ndim, nparams, values=None, constants=None, params=None, scalar=True, update_pos=None)

Bases: object

__init__(code, npos, ndim, nparams, values=None, constants=None, params=None, scalar=True, update_pos=None)

Compiles the matrix multiplication GPU kernel

Parameters

• code – string with CUDA code to evaluate single element of the dense matrix

• npos – number of positions where the field should be evaluated

• ndim – number of dimensions (2 or 3)

• nparams – number of parameters in the cost function

• values – dict of keywords that will be formatted into the code string

• constants – dict of read-only arrays that should be uploaded to the GPU

• params – dict of read/write arrays that should be uploaded to the GPU

• scalar – if True each dimension of the vector field is computed independently

• update_pos – function to update the positions. Gets passed in all the parameter GPU Arrays and the new positions

code(derivative, inverse)

Formats the code with the values embedded

compile_func(derivative, inverse)

Compiles the matrix multiplication kernel on the GPU

Parameters

• derivative – whether to compute the derivative of the field (-1 for no derivative)

• inverse – if True go from vector field derivatives to parameter derivatives

Returns

function that calls the GPU kernel

__call__(parameters, inverse=False, derivative=-1)

Runs matrix multiplication on the GPU

Parameters

• parameters – (nparams, ) array if not inverse else (npos, ndim) array

• inverse – if True go from vector field derivatives to parameter derivatives rather than other way around

• derivative – whether to evaluate derivative of the field (default: no derivative)

Returns

GPUArray with the result

class gyral_structure.cuda.GridConvolve(param_to_field, field_to_param, kernels)

Bases: object

Convolves values in a grid with 4x4x4 convolutions.

__init__(param_to_field, field_to_param, kernels)

Sets up a new grid convolution.

Parameters

• param_to_field – (Np, 4, 4, 4) int array; relates each parameter index to the field in a 4x4x4 neighbourhood; -1 for unused field indices

• field_to_param – (Nf, 4, 4, 4) int array; relates each field index to the parameters in a 4x4x4 neighbourhood

• kernels – (N, 4, 4, 4, 3) array mapping the model parameters to N field parameters

code(inverse=True)

property n_kernels

property n_params

property n_field

__call__(params, inverse=False)

Call self as a function.

class gyral_structure.cuda.QuadraticConvolution(param_to_field, field_to_param, kernel)

Bases: object

Maps the parameters to a single value.

For parameters p_{ijkd} at position (i, j, k) in the 4x4x4 neighbourhood and dimension d, compute sum_{ijkd i’j’k’d’} K_{ijkd i’j’k’d’} a_{ijkd} a_{i’j’k’d’}

__init__(param_to_field, field_to_param, kernel)

Initialize self. See help(type(self)) for accurate signature.

code()

code_derivative_convolve()

code_derivative_sum()

run_gpu(params)

Evaluates the convolution for the given parameters.

The result is stored in _gpu_field.

Parameters

params – numpy array giving the parameters setting the vector potential.

run_gpu_der(params)

Evaluates the derivative of the convolution for the given parameters.

The intermediate result is stored in _gpu_params_allder. The final result is stored in _gpu_params

Parameters

params – numpy array giving the parameters setting the vector potential

per_voxel(params)

Evaluates the convolution for the given parameters.

Parameters

params – numpy array giving the parameters setting the vector potential.

Returns

(nfield, ) array with the result of the convolution for every voxel

evaluate(params)

Evaluates the convolution for the given parameters.

Parameters

params – numpy array giving the parameters setting the vector potential.

Returns

float with the total result of the convolution

derivative(params)

Evaluates the derivative of the quadratic convolution for the given parameters.

Parameters

params – numpy array giving the parameters setting the vector potential.

Returns

numpy array with the same shape of params giving the derivative of the parameters.

gyral_structure.cuda.cached_func

gyral_structure.cuda.cached_source_module

gyral_structure.request: define which voxels/vertices to evaluate

Defines where the vector field (and hence the cost function) will be evaluated.

These are all derived from the FieldRequest. The options are:

• PositionRequest: provide (npos, ndim) array where the field will be evaluated. In most places, where a request is expected a (npos, ndim) array can be provided again, which will be converted into a PositionRequest automatically

• VoxelRequest: provide (npos, ndim) array with voxel centres and scalar or (ndim, ) array with voxel sizes. Will average the field over the voxel by evaluating the field in the center of `nsplit`^3 sub-voxels.

• VertexRequest: provide (npos, ndim, 3) array with triangles. Will average the field over each triangle by evaluating the field in the center of `nsplit`^2 sub-triangles

• DerivativeRequest: provide (npos, ndim) array and decide along which dimension the field derivative will be computed at those points.

• ParamRequest: singular instance of request that can be used to simply return the parameters unaltered rather than computing a vector field

Under the hood to evaluate the field at multiple positions when cost functions are combined a MultRequest is used. The code there ensures that each request is only evaluated once for a given set of parameter even it is used in multiple cost functions.

exception gyral_structure.request.NoCombineError

Bases: Exception

Raise when two requests can not be combined into a single request.

class gyral_structure.request.FieldRequest

Bases: object

Generic object requesting a certain field to be computed.

Can be passed to cost functions (see e.g. cost.VectorCost) to constrain the field solution.

positions = None

get_evaluator(basis, method=())

Gets a function or matrix to compute the requested positions.

Parameters

• basis – Basis function used to compute the field

• method – Computational method used to compute the field

property npos

property ndim

bounding_box(extent=0)

Computes a bounding box around all the positions required by this request.

Parameters

extent – float or (ndim, ) or (2, ndim) array with how much to extend the bounding box beyond the extreme points

Returns

((min(x), max(x)), (min(y), max(y)), (min(z), max(z))) covering the full dataset

radius()

Maximal extent of the request from its center

Returns

scalar or (npos, ) array with the maximal distance from the center.

center()

Center of mass of the request

Returns

(npos, 3) array with the center of mass

shift()

Returns the shift needed so that the positions are centered on zero.

scaling(as_scalar=False)

Returns the extent of the data in each dimension.

Parameters

as_scalar – returns the maximum across all dimensions if True

as_volume(resolution, extent=0)

Defines a volume spanning the complete model area.

Parameters

• resolution – float or (ndim, ) float array with the volume resolution

• extent – float or (ndim, ) or (2, ndim) float array with the distance in mm to go beyond the constraints in each dimension.

Returns

tuple with the array shape and affine transformation

combine(other)

flatten()

dict_to_field(field_dict)

field_to_dict(field)

wrap_qp(qp, Amat)

split()

Splits the Request into positions where the basis function can be evaluated.

Returns

iterable yielding (positions, weights) pairs where the field requested is the weighted sum of the field at all the positions in the iterable

split_cuda()

Returns the code needed to evaluate the positions on the GPU

Returns

tuple with:

• string containing the cuda code in a function named draw_pos

• dictionary with a mapping from variable names to (signature, array) pairs

min_dist(points)

Estimates the minimal distance from a set of points to the requested volume.

Parameters

points – (N, ndim) array for N points in n-dimensional space

Returns

(N, ) array with the minimal distance for every point

class gyral_structure.request.PositionRequest(positions, method=())

Bases: gyral_structure.request.FieldRequest

Computes the field at the provided positions.

Can be passed to cost functions (see e.g. cost.VectorCost) to constrain the field solution.

__init__(positions, method=())

Initialize self. See help(type(self)) for accurate signature.

combine(other)

Combine with another PositionRequest to produce a single PositionRequest that evaluates both

Parameters

other – another PositionRequest

Returns

combined PositionRequest

Raises

NoCombineError if failed to combine

split()

Splits the Request into positions where the basis function can be evaluated.

Returns

iterable yielding (positions, weights) pairs where the field requested is the weighted sum of the field at all the positions in the iterable

split_cuda()

Returns the code needed to evaluate the positions on the GPU

Returns

tuple with:

• string containing the cuda code in a function named draw_pos

• dictionary with a mapping from variable names to (signature, array) pairs

property total_split

class gyral_structure.request.DerivativeRequest(positions, nder, method=(), epsilon=1e-06)

Bases: gyral_structure.request.FieldRequest

Computes the Jacobian numerically at the provided positions.

Can be passed to cost functions (see e.g. cost.VectorCost) to constrain the field solution.

Provides an (npos, ndim, ndim) array where element [:, i, j] contains the derivative of the field component i with respect to the dimension j.

__init__(positions, nder, method=(), epsilon=1e-06)

Initialize self. See help(type(self)) for accurate signature.

static all_derivatives(positions, total_der=1, method=())

Computes all derivatives of a certain power.

Parameters

• positions – (npos, ndim) array where the derivative will be calculated

• total_der – total derivative across all dimensions (defaults to a Jacobian).

• method – overrides the computational method used

combine(other)

radius()

Maximal extent of the request from its center

Returns

scalar or (npos, ) array with the maximal distance from the center.

split_cuda()

Returns the code needed to evaluate the positions on the GPU

Returns

tuple with:

• string containing the cuda code in a function named draw_pos

• dictionary with a mapping from variable names to (signature, array) pairs

split()

Splits the Request into positions where the basis function can be evaluated.

Returns

iterable yielding (positions, weights) pairs where the field requested is the weighted sum of the field at all the positions in the iterable

property total_split

class gyral_structure.request.VoxelRequest(positions, size=1, method=(), nsplit=5)

Bases: gyral_structure.request.FieldRequest

Computes the field averaged over one or more voxels.

Can be passed to cost functions (see e.g. cost.VectorCost) to constrain the field solution.

__init__(positions, size=1, method=(), nsplit=5)

Initialize self. See help(type(self)) for accurate signature.

property nsplit

combine(other)

Combine with another VoxelRequest to produce a single VoxelRequest that evaluates both

Parameters

other – another VoxelRequest

Returns

combined VoxelRequest

Raises

NoCombineError if failed to combine

split()

Splits the Request into positions where the basis function can be evaluated.

Returns

iterable yielding (positions, weights) pairs where the field requested is the weighted sum of the field at all the positions in the iterable

split_cuda()

Returns the code needed to evaluate the positions on the GPU

Returns

tuple with:

• string containing the cuda code in a function named draw_pos

• dictionary with a mapping from variable names to (signature, array) pairs

property total_split

bounding_box(extent=0)

Computes a bounding box around all the positions required by this request.

Parameters

extent – float or (ndim, ) or (2, ndim) array with how much to extend the bounding box beyond the extreme points

Returns

((min(x), max(x)), (min(y), max(y)), (min(z), max(z))) covering the full dataset

radius()

Maximal extent of the request from its center

Returns

scalar or (npos, ) array with the maximal distance from the center.

class gyral_structure.request.VertexRequest(triangles, method=(), nsplit=2)

Bases: gyral_structure.request.FieldRequest

Computes the field averaged over one or more vertices.

Can be passed to cost functions (see e.g. cost.VectorCost) to constrain the field solution.

__init__(triangles, method=(), nsplit=2)

Computes the field averaged over several triangles.

Parameters

• triangles – (npos, ndim, 3) array representing the 3 corner points for npos triangles.

• method – which algorithm to use when evaluating the vertex fields

• nsplit – how often to subdivide the triangles into 6 equal-sized sub-triangles during averaging (so that each field is the average of 6 ** nsplit sub-fields).

property npos

property ndim

property nsplit

property total_split

bounding_box(extent=0)

Computes a bounding box around all the positions required by this request.

Parameters

extent – float or (ndim, ) or (2, ndim) array with how much to extend the bounding box beyond the extreme points

Returns

((min(x), max(x)), (min(y), max(y)), (min(z), max(z))) covering the full dataset

combine(other)

Combine with another VertexRequest to produce a single VertexRequest that evaluates both

Parameters

other – another VertexRequest

Returns

combined VertexRequest

Raises

NoCombineError if failed to combine

size()

Computes the size of the surface elements

split()

Splits the Request into positions where the basis function can be evaluated.

Returns

iterable yielding (positions, weights) pairs where the field requested is the weighted sum of the field at all the positions in the iterable

split_cuda()

Returns the code needed to evaluate the positions on the GPU

Returns

tuple with:

• string containing the cuda code in a function named draw_pos

• dictionary with a mapping from variable names to (signature, array) pairs

radius()

Maximal extent of the request from its center

Returns

scalar or (npos, ) array with the maximal distance from the center.

center()

Center of mass of the request

Returns

(npos, ndim) array with the center of mass

class gyral_structure.request.MultRequest(requests)

Bases: gyral_structure.request.FieldRequest

Container containing multiple requests

The actual requests are stored in the attribute requests

__init__(requests)

Represents one or more requests.

property npos

property ndim

bounding_box(extent=0)

Computes a bounding box around all the positions required by this request.

Parameters

extent – float or (ndim, ) or (2, ndim) array with how much to extend the bounding box beyond the extreme points

Returns

((min(x), max(x)), (min(y), max(y)), (min(z), max(z))) covering the full dataset

flatten()

dict_to_field(field_dict)

field_to_dict(field)

wrap_qp(qp, Amat)

center()

Center of mass of the request

Returns

(npos, 3) array with the center of mass

radius()

Maximal extent of the request from its center

Returns

scalar or (npos, ) array with the maximal distance from the center.

min_dist(points)

Estimates the minimal distance from a set of points to the requested volume.

Parameters

points – (N, ndim) array for N points in n-dimensional space

Returns

(N, ) array with the minimal distance for every point

gyral_structure.request.read_surface(gifti_surf, gifti_mask=None, nsplit=3, method=())

Requests the field to be computed across the surface.

Parameters

• gifti_surf (str) – .surf.gii file describing the cortical surface (typically the white/gray matter boundary)

• gifti_mask – .shape.gii file used to mask the surface (only non-zero vertices will be included; optional)

• nsplit – number of sub-triangles to split each triangle into along each dimension

• method – Algorithm to use when evaluating the basis function at the surface

Returns

Request object asking for the field to be computed across the surface

Return type

VertexRequest

gyral_structure.request.read_volume(nifti_vol, average=True, nsplit=5, in_mm=True, method=())

Requests the field to be computed in all voxels in a volume.

Parameters

• nifti_vol (Union[str, tuple, nibabel.analyze.SpatialImage]) – NIFTI image; any non-zero voxels will be included in the request

• average – set to true to average over the voxel (otherwise the central value is taken).

• nsplit – number of sub-voxels to split each voxel into along each dimension (has no effect if average is False)

• in_mm – returns the positions/sizes in mm rather than voxel coordinates

• method – Algorithm to use when evaluating the basis function at the surface

Returns

Request object asking the field to be computed throughout the ROI (which can be passed on to cost functions)

gyral_structure.request.read_bedpostx(bedpostx_dir, min_volume=0.05, mask=None, average=True, in_mm=True, method=())

Reads in the crossing fibers from bedpostx_dir with a minimum volume fraction min_volume

Parameters

• bedpostx_dir – name of bedpostx directory

• min_volume – minimum volume fraction to include the crossing fiber

• mask – only include voxels with a non-zero value in the mask (has to be in the same space as the diffusion data)

• average – set to true to average over the voxel (otherwise the central value is taken).

• in_mm – if True request positions in mm rather than voxel space

• method – Algorithm to use when evaluating the basis function at the surface

Returns

request to compute the field and (N vectors, N positions, N dimensions) array of target dyads

gyral_structure.streamline: streamline tractography through vector field

Evaluate the streamlines running through a given field.

gyral_structure.streamline.surface_startpoints(triangles, field_density)

Initializes streamlines at the surface based on the field density.

Parameters

• triangles – (npos, ndim, 3) array representing the 3 corner points for npos triangles.

• field_density – (npos, ) array with the total number of streamlines crossing the surface

Returns

(N, ndim) array of starting points

gyral_structure.streamline.volume_startpoints(charge_pos, charge_size, nstream)

Initializes streamlines at a charge distribution.

Streamlines are generated randomly from a uniform distribution in the sphere

Parameters

• charge_pos – centers of the charge spheres

• charge_size – size of the charge spheres

• nstream – float array; number of streamlines per charge sphere

class gyral_structure.streamline.Tracker(orient_func, stop_mask, affine_mask, step_size=0.1, maxstep=30000, stop_surfaces=())

Bases: object

__init__(orient_func, stop_mask, affine_mask, step_size=0.1, maxstep=30000, stop_surfaces=())

Creates a new tracking algorithm using Runga Kutta RK4

Parameters

• orient_func – (npos, 3) array -> orientation; gives the fiber orientation at any point in space

• stop_mask – Nifti1Image; stop the streamline when it leaves the mask

• stop_affine – voxel -> mm transformation for the stop mask

• stop_conditions – List[(npos, 3) array, (npos, 3) array -> (npos, ) bool array]; tests whether a streamline should stop given its previous and current position

track(start, nstore=10, inverse=False)

update_pos(current_pos, inverse=False)

Runge-Kutta RK4 update of position

gyral_structure.streamline.to_vtk(streamlines, color=(1.0, 0.0, 0.0), linewidth=1, renderer=None, window=None, interact=False)

Convert streamlines to a vtkPolyData dataset for rendering in VTK

The streamlines produced by Tracker.track can be used for plotting in VTK

Parameters

• streamlines – list of individual streamlines (given as (N, 3) arrays)

• color –

  RGB color to show the streamlines; Could be one of:

  • (3, ) array to set the same color to all lines

  • (nlines, 3) array to set a different color to every line

  • (npoints, 3) to set a different color to each line segment

• renderer – vtkRenderer to which the lines will be added (default: new one is created

• window – vtkRenderWindow to which the renderer will be added (default: new one is created)

• interact – if True makes the VTK window interactive

Returns

VTK objects containing the streamlines

Return type

tuple[vtkPolyData, vtkRenderer, vtkRenderWindow]

gyral_structure.utils: several utilities

gyral_structure.utils.get_filetype(filename)

Determines whether the filename is zipped based solely on the filename extension

Returns

File object used to open the filename

gyral_structure.utils.dyad_to_mm(dyads, affine)

Converts the dyads from the FSL convention to mm space

Can be inverted using dyad_from_mm().

Parameters

• dyads – (…, 3) array with the dyads as loaded from e.g. a dti_V1 file

• affine – (4, 4) or (3, 3) array with the affine transformation from voxel to mm space (in the latter case excluding translation)

Returns

(…, 3) array with the normalised dyads in mm space

gyral_structure.utils.dyad_from_mm(dyads, affine)

Converts the dyads from mm space to the FSL convention

Can be inverted using dyad_to_mm().

Parameters

• dyads – (…, 3) array with the dyads in mm space

• affine – (4, 4) or (3, 3) array with the affine transformation from voxel to mm space (in the latter case excluding translation)

Returns

(…, 3) array with the normalised dyads in the FSL convention