gyral_structure.basis: defining the basis functions

Defines a set of basis functions that map from input parameters to a vector field at the requested positions

The locations where the field will be evaluated will be defined by the requests (see request).

Most of these basis functions produce divergence-free vector fields (i.e., no terminating streamlines), although some (e.g. ChargeDistribution) are explicitly designed to model the termination of streamlines in selected gray matter regions.

Multiple basis functions can be added together using SumBase.

New basis functions should be derived from BasisFunc.

Concrete implementations of basis functions (in order of decreasing usefulness) are:

• RadialBasis: set of matrix-valued, radial basis function with limited extent

• ChargeDistribution: vector field generated for charges generating streamline endpoints in the gray matter

• Gradient: A smooth distribution of the vector field potential

• Fourier: Periodic vector field potential

• Polynomial: vector field potential given by polynomial

Excellent results can usually be obtained by using ChargeDistribution to generate a reasonable initialisation of the field, which can be further updated using one or more RadialBasis.

gyral_structure.basis.core module

Defines the two central classes defining the basis function, that maps parameters to vector fields

• BasisFunc should be used as a superclass for any parameter to field mapper

• SumBase can be used to create a field by summing multiple other basis functions

In addition a helper function wrap_numba() is defined, which is meant to assist in writing implementations of the basis functions using numba

gyral_structure.basis.core.wrap_numba(npos, ndim, nparams, scalar=True, sparse=None)

Vectorizes a numba function that computes a single element of the matrix mapping parameters to fields

The resulting function can be used to define a basis function

Parameters

• npos – number of positions where the field will be calculated

• ndim – dimensionality of the space (i.e. 2- or 3-dimensional)

• nparams – number of parameters

• scalar – If True the numba function only computes a scalar output rather than the full vector

• sparse – Optionally two equal-length lists of integers indexing which positions and parameters affect each other

Returns

wrapper function that converts an input function with signature (position, idx_param, derivative=-1) -> float to an output function with signature (positions, parameters, derivative=-1, inverse=False) -> field

class gyral_structure.basis.core.BasisFunc(nparams)

Bases: object

Defines a basis function mapping 2/3D positions to a 2/3D vector field.

Call this object with an (npos, ndim) array of positions to get a function which computes the field from parameters or the other way around

The user of any subclass will typically call one of two methods:

• get_evaluator() which returns a RequestEvaluator that maps parameters to the field at a fixed position (used during field optimisation)

• param_evaluator() which returns a function that maps positions to the field for a fixed set of parameters (used during streamline tractography).

Subclasses should override three methods:

• get_partial_mat(): returns the matrix to compute the potential or vector field at given positions

• get_func(): return functions that compute the potential or vector field at given positions

• validate_pos(): returns True for any positions where the basis function has been defined (optional)

What these functions return is defined by setting the following attributes:

• potential

• scalar

• compute

potential = True

If True the matrix/function is expected to map to the vector potential rather than vector field (ensures that the result is divergence-free)

scalar = True

If True the matrix/function is expected to return only a single dimension o the vector field or potential rather than the full vector field or potential

ndim = 3

dimensionality of the space in which the basis function is defined

compute = ()

sequence of possible algorithms that can be used in the evaluation (e.g., is there a gpu or a matrix implementation?)

__init__(nparams)

Creates a new Basis function.

Parameters

nparams – number of parameters that will be passed to the basis function

property scalar_nparams

Number of parameters used to define a single dimension of the output vector field

get_partial_mat(positions, derivative=-1)

Computes the matrix used to map parameters to the vector field or potential at the given positions

Should be overwritten by sub-classes. What it exactly calculates depends on the potential and scalar flags:

• if the potential flag is set (default) then the matrix should compute the vector potential rather than fiber orientations

• if the scalar flag is set return one value per position rather than ndim

Parameters

• positions – (npos, 3) array with the positions of interest

• derivative – which spatial derivative to compute (default: no derivative)

Returns

(npos, nparams // 3) array mapping the parameters to the positions

stack_matrices(positions)

Returns the (npos, nparams) matrix that maps the parameters to the n-dimensional field.

Combines the sub-matrices returned by get_partial_mat on the given positions.

Parameters

positions – (npos, ndim) array with the locations where the field should be evaluated

Returns

(npos, nparams) matrix mapping the parameters to the vector field

get_full_mat(req)

Returns the matrix used to __call__ the given request.

Parameters

req (request.FieldRequest) – Defines where the field should be evaluated

Returns

(npos, nparams) matrix mapping the parameters to the vector field

get_func(positions, method)

Creates a function that evaluates the requested field at the positions.

Should be overwritten by sub-classes. What it exactly calculates depends on the potential and scalar flags: - if the potential flag is set (default) then the matrix should compute the vector potential rather than fiber orientations - if the scalar flag is set return one value per position rather than ndim

Parameters

• positions – (npos, ndim) array defining where the result should be computed

• method – which algorithm to use in the computation

Returns

function that maps (input, inverse, derivative) -> output (field derivative to parameter derivative if inverse else param to field)

Raises

NotImplementedError if the basis function is not defined for a specific method

get_evaluator(req, method=None)

Computes the function used to __call__ the request for a given method.

raises NotImplementedError if the method has not been implemented for this basis function and request.

Parameters

• req (request.FieldRequest) – defines where the field should be evaluated

• method (algorithm.Algorithm or string) – the preferred algorithm to use when evaluating the field at the requested positions

Returns

an RequestEvaluator that has been set up to quickly evaluate the field at the requested positions given a vector field.

param_evaluator(parameters, method=None, extent=0, nsim=4096)

Returns a function that evaluates the field at a given location

Used for streamline evaluation.

Parameters

• parameters – set of parameters for which the field should be evaluated at different positions

• method – algorithm to use when computing field

• extent – maximum extent of request to accept (keep at zero for tractography)

• nsim – maximum number of positions to evaluate simultaneously

Returns

function to map the positions to the field

validate(req)

Checks for which voxels/vertices/positions in the request the basis function is defined.

Parameters

req (request.FieldRequest) – defines where the basis function should be evaluated

Returns

boolean array which is True, wherever the request can be safely evaluated

validate_pos(positions)

Whether the summed basis function can be computed at the provided positions

Parameters

positions – (npos, 3) array with the positions of interest

Returns

(npos, ) boolean array, which is True if the positions are within range of all basis functions

save(filename, solution=())

Saves the basis function with the solution to an HDF5 file

Parameters

• filename – HDF5 (.h5) file to store the basis function & solution in

• solution – array with the best-fit parameters

precompute_evaluator(extent=0)

Prepares the evaluator generator to quickly generate results for different requests

Called before doing tractography

Parameters

extent – maximum size of the request

to_hdf5(group: h5py._hl.group.Group)

Stores the basis function in the provided HDF5 group

Parameters

group – HDF5 group, where the basis function will be stored

classmethod from_hdf5(group: h5py._hl.group.Group)

Reads the basis function from the provided HDF5 group

Parameters

group – HDF5 group, where the basis function is stored

Returns

newly read basis function

class gyral_structure.basis.core.SumBase(initial_list=None)

Bases: collections.UserList

Combination of multiple basis functions, which are added together to create a single net vector field

After creation the user will typically call one of two methods:

• get_evaluator() which returns a MultRequestEvaluator that maps parameters to the net field at a fixed position (used during field optimisation)

• param_evaluator() which returns a function that maps positions to the net field for a fixed set of parameters (used during streamline tractography).

If not all basis functions should be updated during the optimisation, those can be fixed to a set of parameters using fix() or fix_all(). They can later be released using release().

__init__(initial_list=None)

Combines the fields from one or more Basis functions.

Parameters

initial_list – list of sub-classes of BasisFunc (or SumBase)

property compute

A sequence of algorithms that can be used to define all the basis functions comprising the net field

property nparams

Total number of parameters required to supply all the basis functions

property ndim

Dimensionality of the space in which the basis functions have been defined

fix(index, params)

Fixes a basis function with the given parameters

Fixed basis functions will only compute the resulting vector field once (for the given parameters) and then consistently return that output

Parameters

• index – integer index indicating which basis function will be fixed

• params – array with the parameters used to fix the vector field of the given basis function

fix_all(params)

Fixes all basis function with the given parameters

Fixed basis functions will only compute the resulting vector field once (for the given parameters) and then consistently return that output

Parameters

params – (nparams, ) array with the parameters used to fix the vector field of the given basis function

release(index)

Removes the fix from a given basis function

Fixed basis functions will consistently return the same vector field without getting any parameters

release_all()

Removes the fix from all the basis functions

Fixed basis functions will consistently return the same vector field without getting any parameters

param_evaluator(parameters, method=None, extent=0, nsim=4096)

Returns a function that evaluates the basis function for arbitrary positions given a fixed parameter array

Used to __call__ streamlines

Parameters

• parameters – (nparams, ) array defining the parameters for which the field will be evaluated

• method – Algorithm used to __call__ the basis functions

• extent – maximum extent of request to accept (keep at zero for tractography)

• nsim – maximum number of positions to evaluate simultaneously

Returns

function that maps positions to vector field

validate_pos(positions)

Whether the summed basis function can be computed at the provided positions

Parameters

positions – (npos, 3) array with the positions of interest

Returns

(npos, ) boolean array, which is True if the positions are within range of all basis functions

save(filename, solution=())

Saves the basis functions with the solution to an HDF5 file

Parameters

• filename – HDF5 (.h5) file to store the basis function & solution in

• solution – array with the best-fit parameters

get_evaluator(req, method=None)

Computes the function used to __call__ the request for a given method.

raises NotImplementedError if the method has not been implemented for this basis function and request.

Parameters

• req (request.FieldRequest) – defines where the field should be evaluated

• method (Algorithm) – method used to __call__ the request

to_hdf5(group: h5py._hl.group.Group)

Stores the basis functions in the provided HDF5 group

Parameters

group – HDF5 group, where the basis function will be stored

classmethod from_hdf5(group: h5py._hl.group.Group)

Reads the basis functions from the provided HDF5 group

Parameters

group – HDF5 group, where the basis function is stored

Returns

newly read basis function

gyral_structure.basis.core.read(filename)

Loads the basis function and best-fit parameters from a pickle file.

Parameters

filename – file to load the basis function & solution from

gyral_structure.basis.core.read_pickle(filename)

Loads the basis function and best-fit parameters from a pickle file.

Pickle format is now outdated, use HDF5 instead

Parameters

filename – file to load the basis function & solution from

gyral_structure.basis.core.get_basisfunc(name)

Gets the basis function from the class name

Parameters

name – name of a subclass of Basisfunc or SumBase

Returns

class with the name

Raise

KeyError if subclass does not exist (or has not been imported yet)

gyral_structure.basis.evaluator module

This module defines the mapping from parameters to vector field given a basis function (from basis) and a request defining where the field should be evaluated (from request).

For a basis function and request a new evaluator can be obtained using get_single_evaluator(). The result will be a subclass of RequestEvaluator.

gyral_structure.basis.evaluator.get_single_evaluator(basis, req)

Main function to call to get an evaluator of a req

The appropriate evaluator is chosen based on the first match in the following rules:

• ParamRequest -> IdentityEvaluator

• Multiple requests or basis functions -> MultEvaluator

After this we iterate through the available methods to find a valid algorithm (in order we check the methods provided when calling get_evaluator, the methods provided to the req, the methods provided to the basis function, the possible algorithms for the basis function):

• first valid algorithm is matrix or matrix_cuda -> MatRequestEvaluator

• otherwise -> FuncRequestEvaluator

Parameters

• basis (core.BasisFunc) – basis function mapping parameters to field

• req (req.FieldRequest) – defines the positions/voxels/vertices where the field should be evaluated

Returns

Callable object that converts parameters into the field at the requested position (or the opposite for field derivatives)

class gyral_structure.basis.evaluator.RequestEvaluator(basis, req, method)

Bases: object

Evaluates the basis function at the requested positions

The user interaction with this object will usually consist of calling it:

• evaluator(params) maps from the parameters to the vector field

• evaluator(field_derivatives, inverse=True) maps from the field derivatives to the parameter derivatives

For efficiency this will call under the hood:

• evaluate_results(): initializes all the calculations

• combine_results(): combines the results of the calculations

• clean_results(): cleans the previous results

__init__(basis, req, method)

Creates a new basis function to __call__() the request

Parameters

• basis (core.BasisFunc) – basis function mapping parameters to field

• req (request.FieldRequest) – defines the positions/voxels/vertices where the field should be evaluated

• method (Algorithm) – defines the algorithm used to __call__ the basis function

__call__(params, inverse=False)

Converts a set of parameters to the vector field at the requested positions (or inverse of that)

Parameters

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

• inverse – set to true to convert a vector field derivative into parameter derivatives

Returns

(nparams, ) array if inverse else (npos, ndim) array

evaluate_results(params, inverse=False)

Runs the basis function for all positions in the request without accessing the results

The results of this can be obtained from combined_results after running evaluate_results

Parameters

• params – (npos, dim) array if inverse else (nparams, ) array

• inverse – set to true to convert a vector field derivative into parameter derivatives

combine_results(inverse=False)

Retrieves the vector field at the requested positions as computed by evaluate_results

evaluate_results should be called with the same parameters before running this command

Parameters

inverse – set to true to convert a vector field derivative into parameter derivatives

clean_results()

Cleans previously evaluated results

Can be used to free the memory used up by the previous results

property scalar

Maps to self.basis.scalar

property potential

Maps to self.basis.potential

property ndim

Maps to self.basis.ndim

property npos

Maps to self.request.npos

property nparams

Maps to self.basis.nparams

property use_mat

If True pre-generate a matrix rather than evaluating it on the fly

property use_cuda

If True evaluates the field on the GPU rather than CPU

update_pos(new_request)

class gyral_structure.basis.evaluator.IdentityEvaluator(basis, req, method=None)

Bases: gyral_structure.basis.evaluator.RequestEvaluator

Returns the parameters irrespective of the basis function

use_cuda = False

__init__(basis, req, method=None)

Creates a new basis function to __call__() the request

Parameters

• basis (core.BasisFunc) – basis function mapping parameters to field

• req (request.FieldRequest) – defines the positions/voxels/vertices where the field should be evaluated

• method (Algorithm) – defines the algorithm used to __call__ the basis function

evaluate_results(params, inverse=False)

Runs the basis function for all positions in the request without accessing the results

The results of this can be obtained from combined_results after running evaluate_results

Parameters

• params – (npos, dim) array if inverse else (nparams, ) array

• inverse – set to true to convert a vector field derivative into parameter derivatives

combine_results(inverse=False)

Retrieves the vector field at the requested positions as computed by evaluate_results

evaluate_results should be called with the same parameters before running this command

Parameters

inverse – set to true to convert a vector field derivative into parameter derivatives

clean_results()

Cleans previously evaluated results

Can be used to free the memory used up by the previous results

__call__(params, inverse=False)

Returns the input parameters

update_pos(new_request)

class gyral_structure.basis.evaluator.FuncRequestEvaluator(basis, req, method)

Bases: gyral_structure.basis.evaluator.RequestEvaluator

Evaluates the basis function using on-the-fly evaluation of the matrix on the CPU or GPU

The user interaction with this object will usually consist of calling it:

• evaluator(params) maps from the parameters to the vector field

• evaluator(field_derivatives, inverse=True) maps from the field derivatives to the parameter derivatives

For efficiency this will call under the hood:

• evaluate_results(): initializes all the calculations

• combine_results(): combines the results of the calculations

• clean_results(): cleans the previous results

__init__(basis, req, method)

Creates a new basis function to __call__() the request

Parameters

• basis (core.BasisFunc) – basis function mapping parameters to field

• req (request.FieldRequest) – defines the positions/voxels/vertices where the field should be evaluated

• method (Algorithm) – defines the algorithm used to __call__ the basis function

func(params, inverse=False, derivative=-1)

Calls the basis function for all positions in the request

evaluate_results(params, inverse=False)

Runs the basis function for all positions in the request without accessing the results

The results of this can be obtained from combined_results after running evaluate_results

Parameters

• params – (npos, dim) array if inverse else (nparams, ) array

• inverse – set to true to convert a vector field derivative into parameter derivatives

combine_results(inverse=False)

Retrieves the vector field at the requested positions as computed by evaluate_results

evaluate_results should be called with the same parameters before running this command

Parameters

inverse – set to true to convert a vector field derivative into parameter derivatives

clean_results()

Cleans previously evaluated results

Can be used to free the memory used up by the previous results

update_pos(new_request)

class gyral_structure.basis.evaluator.MatRequestEvaluator(basis, req, method)

Bases: gyral_structure.basis.evaluator.RequestEvaluator

Evaluates the basis function with a pre-defined matrix on CPU or GPU

The user interaction with this object will usually consist of calling it:

• evaluator(params) maps from the parameters to the vector field

• evaluator(field_derivatives, inverse=True) maps from the field derivatives to the parameter derivatives

For efficiency this will call under the hood:

• evaluate_results(): initializes all the calculations

• combine_results(): combines the results of the calculations

• clean_results(): cleans the previous results

__init__(basis, req, method)

Creates a new basis function to __call__() the request

Parameters

• basis (core.BasisFunc) – basis function mapping parameters to field

• req (request.FieldRequest) – defines the positions/voxels/vertices where the field should be evaluated

• method (Algorithm) – defines the algorithm used to __call__ the basis function

evaluate_results(params, inverse=False)

Runs the matrix multiplication and temporarily stores the result

The results of this can be obtained from combined_results after running evaluate_results

Parameters

• params – (npos, dim) array if inverse else (nparams, ) array

• inverse – set to true to convert a vector field derivative into parameter derivatives

combine_results(inverse=False)

Retrieves the vector field at the requested positions as computed by evaluate_results

evaluate_results should be called with the same parameters before running this command

Parameters

• params – (npos, dim) array if inverse else (nparams, ) array

• inverse – set to true to convert a vector field derivative into parameter derivatives

clean_results()

Cleans previously evaluated results

Can be used to free the memory used up by the previous results

wrap_qp(qp)

update_pos(new_request)

class gyral_structure.basis.evaluator.MultEvaluator(sum_basis, full_request)

Bases: object

Evaluates multiple basis functions/requests

The user interaction with this object will usually consist of calling it:

• evaluator(params) maps from the parameters to the vector field

• evaluator(field_derivatives, inverse=True) maps from the field derivatives to the parameter derivatives

For efficiency this will call under the hood:

• evaluate_results(): initializes all the calculations

• combine_results(): combines the results of the calculations

• clean_results(): cleans the previous results

__init__(sum_basis, full_request)

Creates a new basis function/request combination

Parameters

• sum_basis – basis function (optionally sum of multiple)

• full_request – request (optionally combination of multiple)

fixed_field = None

Total field from all the basis functions, whose parameters have been fixed

__call__(params, inverse=False)

Converts a set of parameters to the vector field at the requested positions (or inverse of that)

Parameters

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

• inverse – set to true to convert a vector field derivative into parameter derivatives

Returns

(nparams, ) array if inverse else (npos, ndim) array

static loop_basis(basis_func, current=0)

Loops over all the basis functions to find those that need to be evaluated and the fixed field

Parameters

• basis_func – SumBasis collection of basis function or any individual basis function

• current – keep this at zero to get an accurate number of expected parameters

Returns

Tuple with:

• list of (basis function, slice of parameters that should enter basis function)

• integer with total number of expected parameters

• (npos, ndim) array with the field for the fixed basis functions

get_evaluator(basis, request)

evaluate_results(params, inverse=False)

Call evaluate_results on all evaluators with the correct parameters

combine_results(inverse=True)

clean_results()

mat()

wrap_qp(qp)

property use_cuda

property use_mat

update_pos(new_positions)

gyral_structure.basis.grid module

This module is no longer actively maintained

class gyral_structure.basis.grid.Grid(mask, affine=None)

Bases: gyral_structure.basis.core.BasisFunc

Basis function defined on a grid.

Every grid element should only affect the field in the direct neighbourhood.

ndim = 3

params_per_grid = None

__init__(mask, affine=None)

Creates a grid at least covering the voxels in the mask.

The field will be computable in every voxel centered on the locations of the non-zero values in mask

Parameters

• mask – True where the field should be computable

• affine – conversion from voxel to mm space (default: no transformation)

property nparams

Total number of parameters describing the vector potential

to_indices(positions, return_rel_pos=False)

Returns the voxel indices that the positions are on the self.idx_field grid

Parameters

• positions – (npos, ndim) array of positions of interest in mm

• return_rel_pos – if True, returns the position relative to the voxel center (between -0.5 and 0.5) as well

Returns

(npos, ndim) int array of voxel indices. Also returns an (npos, ndim) float array of offsets from the voxel center.

validate_pos(positions)

Tests whether the field is defined at the positions

Parameters

positions – (npos, ndim) array of positions of interest in mm

Returns

(npos, ) bool array, which is True if the position is on the mask

contributors(positions)

Returns the indices of the parameter indices contributing to defining the field at positions

Parameters

positions – (npos, ndim) float array of positions of interest in mm

Returns

(npos, 2, 2, 2) int array defining the parameter indices describing the potential in the (2x2x2) neighbourhood around the positions

get_full_mat(positions)

Returns the (npos, nparams) matrix that maps the parameters to the 3-dimensional field in mm space.

Combines the sub-matrices returned by get_partial_mat on the positions.

The function in BasisFunc is overriden so that the field affine transformation can be applied to map the derivative in voxel space to the derivative in mm space

Parameters

positions – (npos, ndim) positions where the field should be evaluated

get_full_func(positions, method)

Combines the evaulations provided by get_func to create a function that maps the parameters to the full field at positions

The function in BasisFunc is overriden so that the field affine transformation can be applied to map the derivative in voxel space to the derivative in mm space

Parameters

• positions – (npos, ndim) array defining where the result should be computed

• method – which algorithm to use in the computation

Returns

function that maps (input, inverse) -> output (derfield to derparam if inverse else param to field)

apply_affine(orientations, inverse=False)

Converts orientations in voxel-space to mm-space (by default) or from mm-space to voxel-space (if inverse is True)

NOTE: fiber orientations are expected to be in mm-space not voxel-space

Parameters

• orientations – (norient, ndim) array defining the orientations in voxel space (or mm space if inverse is True)

• inverse – if set convert from mm to voxel space rather than the other way around

Returns

(norient, ndim) array defining the orientations in mm space (or voxel space if inverse is True)

class gyral_structure.basis.grid.GridFuncRequestEvaluator(basis, req, method)

Bases: gyral_structure.basis.evaluator.FuncRequestEvaluator

evaluate_results(params, inverse=False)

Runs the basis function for all positions in the request without accessing the results

The results of this can be obtained from combined_results after running evaluate_results

Parameters

• params – (npos, dim) array if inverse else (nparams, ) array

• inverse – set to true to convert a vector field derivative into parameter derivatives

combine_results(inverse=False)

Retrieves the vector field at the requested positions as computed by evaluate_results

evaluate_results should be called with the same parameters before running this command

Parameters

inverse – set to true to convert a vector field derivative into parameter derivatives

gyral_structure.basis.radial module

Defines a basis function based on compactly supported radial basis function.

Two options are provided to define the shape of the radial basis function:

• wendland()

• buhmann()

The function hcp_packing() can be used to define an hexagonal packing configuration for the RBF centroids.

A choice of centroids and RBF shape are needed to create the basis function RadialBasis.

gyral_structure.basis.radial.hcp_packing(bounding_box, distance=1)

Proposes an hexagonal packing for the centroids of a radial basis function.

Parameters

• bounding_box – ((min(x), max(x)), (max(y), max(y)), (max(z), max(z))) covering the full data-set

• distance – minimum distance between any two centroids

Returns

(npos, 3) array with the proposed centroids

Return type

sp.ndarray

class gyral_structure.basis.radial.RadialBasis(mapping, centroids, size=1.0)

Bases: gyral_structure.basis.core.BasisFunc

Basis function for compactly supported any radial basis functions

The radial basis function used is the one provided by mapping.

potential = False

RBFs directly calculate the full vector field, not the vector field potential

scalar = False

The mapping from parameters to vector field is not independent across dimensions

compute = ('cuda', 'matrix_cuda', 'matrix', 'numba')

RBFs have been implemented for the GPU and CPU and both matrix and on-the-fly format

__init__(mapping, centroids, size=1.0)

Creates a new set of basis functions.

Parameters

• mapping (CompactRBFSympy) – callable mapping the offsets from the centroids to fields (from buhmann() or wendland())

• centroids – (nparams, ndim) array of the centroids of the radial basis functions (e.g., from hcp_packing())

• size – (nparams, ) array or scalar; total extent of the radial basis functions

get_full_mat(req)

Returns the matrix used to __call__ the given request.

Parameters

req (request.FieldRequest) – Defines where the field should be evaluated

Returns

(npos * ndim, nparams) matrix mapping the parameters to the vector field

get_partial_mat(positions, derivative=-1)

Computes the matrix used to map parameters to the vector field or potential at the given positions

Parameters

• positions – (npos, 3) array with the positions of interest (shift and scaling have already been applied)

• derivative – which spatial derivative to compute (default: no derivative)

Returns

(npos, nparams) array mapping the parameters to the positions

get_func(positions, method)

Creates a function that evaluates the requested field at the positions.

Parameters

• positions – (npos, ndim) array defining where the result should be computed

• method – which algorithm to use in the computation

Returns

function that maps (input, inverse, derivative) -> output (field derivative to parameter derivative if inverse else param to field)

Raises

NotImplementedError if the basis function is not defined for a specific method

property nparams

Number of parameters expected

property ndim

Dimensionality of the space in which the RBF is embedded

within_range(req, return_compressed=False)

Return for every centroid which elements of the request are within range.

For speed during tractography intersections can be precomputed on a grid using precompute_grid

Parameters

• req (request.FieldRequest) – describes where the basis function should be evaluated

• return_compressed – return the request indices in a compressed format

Returns

tuple with the request and centroid indices in compressed format

bounding_box(extent=0)

Computes a bounding box around all the centroids where requests might be affected

Parameters

extent – float with the maximum size of the request

Returns

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

precompute_evaluator(extent=0, resolution=None)

Precomputes which dipoles intersect with a grid

The result is stored internally and will be used when calling within_range().

Parameters

• extent – maximum size of the request

• resolution – resolution of the grid on which the intersections will be precomputed

to_hdf5(group: h5py._hl.group.Group)

Stores the radial basis function in the provided HDF5 group

Parameters

group – HDF5 group, where the basis function will be stored

classmethod from_hdf5(group: h5py._hl.group.Group)

Reads the radial basis function from the provided HDF5 group

Parameters

group – HDF5 group, where the basis function is stored

Returns

newly read radial basis function

gyral_structure.basis.radial.wendland(smoothness=2)

Compactly-supported radial basis functions from Wendland (1995)

Parameters

smoothness – Sets how smooth the radial basis function is: 0: (1 - r)+^2 1: (1 - r)+^4 (4 r + 1) 2: (1 - r)+^6 (35 r^2 + 18 r + 3) 3: (1 - r)+^8 (32 r^3 + 25 r^3 + 8r + 1)

Returns

object that evaluates the radial basis function

Return type

CompactRBFSympy

gyral_structure.basis.radial.buhmann(n=3)

Compactly-supported radial basis functions from Buhmann (2001)

Parameters

n – Chooses from the RBFs explicitly stated by Buhmann: 1: 1/18 - r^2 + 4/9 r^3 + 1/2 r^4 - 4/3 r^3 log(r) (twice differentiable) 2: 2 r^4 log r - 7/2 r^4 + 16/3 r^3 - 2 r^2 + 1/6 (twice differentiable) 3: 112 / 45 r^(9/2) + 16 / 3 r^(7/2) - 7 r^4 - 14 / 15 r^2 + 1/9 (thrice differentiable)

Returns

object that evaluates the radial basis function

Return type

CompactRBFSympy

class gyral_structure.basis.radial.CompactRBFSympy(expr)

Bases: object

Sympy representations of compactly supported radial basis function

__init__(expr)

Creates a new RBF from the given sympy expression

Parameters

expr – sympy expression defining density as a function of ‘r’

__call__(offset, radius=None)

Evaluates the basis function

Parameters

• offset – (N, 3) array with the offset between the centroids and the requested positions

• radius – (N, ) with the size of the offset. Will be computed if not provided

Returns

(N, ndim, ndim) array with the effect size of any given parameter on the field

get_func(positions, centroids, size, method)

Creates a function that evaluates the requested field at the positions.

Parameters

• positions – (npos, ndim) array defining where the result should be computed

• centroids – (ncen, ndim) array defining where the radial basis functions are located

• size – (ncen, ) array with the extent of the radial basis functions

• method – which algorithm to use in the computation

Returns

function that maps (input, inverse, derivative) -> output (field derivative to parameter derivative if inverse else param to field)

Raises

NotImplementedError if the basis function is not defined for a specific method

gyral_structure.basis.radial.RadialGPUArrays = <gyral_structure.basis.radial._RadialGPUArrays object>

Singleton object with the Radial GPU arrays used in the evaluator

class gyral_structure.basis.radial.RadialBasisCudaEvaluator(basis, request, method=None)

Bases: gyral_structure.basis.evaluator.RequestEvaluator

Specialised evaluator to compute the radial basis functions on the GPU

global_cuda_params = {}

use_cuda = True

use_mat = False

__init__(basis, request, method=None)

CUDA evaluator for radial basis functions.

Maps parameters to vector field or the inverse

Parameters

• basis (RadialBasis) – Radial basis function

• request (request.FieldRequest) – Any single request

• method – overrides the Computational method used

update_indices(request=None, only_forward=False)

evaluate_results(params, inverse=False)

Runs the basis function for all positions in the request without accessing the results

The results of this can be obtained from combined_results after running evaluate_results

Parameters

• params – (npos, dim) array if inverse else (nparams, ) array

• inverse – set to true to convert a vector field derivative into parameter derivatives

combine_results(inverse=False)

Retrieves the vector field at the requested positions as computed by evaluate_results

evaluate_results should be called with the same parameters before running this command

Parameters

• params – (npos, dim) array if inverse else (nparams, ) array

• inverse – set to true to convert a vector field derivative into parameter derivatives

clean_results()

Cleans previously evaluated results

Can be used to free the memory used up by the previous results

update_pos(new_request)

gyral_structure.basis.tricubic module

This module is no longer actively maintained

Tricubic interpolation based on the vector potential and jacobians at the grid points.

class gyral_structure.basis.tricubic.Tricubic(mask, affine=None, alldiff=False)

Bases: gyral_structure.basis.grid.Grid

Tricubic interpolation based on the vector potential and the jacobian defined at the voxel corners.

scalar = True

potential = True

property params_per_grid

__init__(mask, affine=None, alldiff=False)

Creates a grid at least covering the voxels in the mask.

The field will be computable in every voxel centered on the locations of the non-zero values in mask

Parameters

• mask – True where the field should be computable

• affine – conversion from voxel to mm space (default: no transformation)

• alldiff –

get_partial_mat(positions, derivative=-1)

Returns the sparse matrix that computes the derivative of a scalar potential at the given positions

Note that the resulting derivatives are in voxel-space, not mm-space (get_full_mat converts them to mm-space).

Parameters

• positions – (npos, ndim) float array of positions of interest in mm

• derivative – -1 to compute the potential, (0, 1, 2) to compute the derivative in respectively (x, y, z)

gyral_structure.basis.tricubic.poly_to_pos_scalar(rel_pos, nder=(0, 0, 0))

Maps the polynomial parameters to the derivative of a scalar field.

Parameters

• rel_pos – (…, 3) array with the relative position wrt voxel center (between -0.5 and 0.5)

• nder – power of the derivative in each spatial dimension

Returns

(…, 64) array mapping the 32 polynomial parameters to the scalar field derivatives

gyral_structure.basis.tricubic.grid_to_pos_scalar(rel_pos, nder=(0, 0, 0))

Maps the scalar/derivative of the scalar potential from the voxel corners to any points within the voxel

Parameters

• rel_pos – (…, 3) array with the relative position wrt voxel center (between -0.5 and 0.5)

• nder – power of the derivative in each spatial dimension

Returns

(…, 2, 2, 2, 8) array with - first 3 dimensions determining which corner of the voxel is considered - 4th dimension determines the derivative (0: no derivative, 1: x, 2: y, 3: z)

gyral_structure.basis.trilinear module

This module is no longer actively maintained

Implements a basis function that defines the potential on a trilinear grid.

class gyral_structure.basis.trilinear.TriLinear(mask, affine=None)

Bases: gyral_structure.basis.grid.Grid

Defines the potential as a trilinear interpolation from a regular parameter grid.

params_per_grid = 3

compute = ('numba', 'matrix', 'matrix_cuda')

get_partial_mat(positions, derivative=-1)

Returns the sparse matrix that computes the derivative of a scalar potential at the given positions

Note that the resulting derivatives are in voxel-space, not mm-space (get_full_mat converts them to mm-space).

Parameters

• positions – (npos, ndim) float array of positions of interest in mm

• derivative – -1 to compute the potential, (0, 1, 2) to compute the derivative in respectively (x, y, z)

get_func(positions, method)

Creates a function that evaluates the derivative of a scalar potential at the positions.

Note that the resulting derivatives are in voxel-space, not mm-space (get_full_func converts them to mm-space).

Parameters

• positions – (npos, ndim) array defining where the result should be computed

• method – which algorithm to use in the computation

Returns

function that maps (input, inverse, derivative) -> output (derfield to derparam if inverse else param to field)

gyral_structure.basis.trilinear.poly_to_pos_scalar(rel_pos, derivative=(0, 0, 0))

Given a position relative to the voxel center computes the mapping from the polynomial parameter to the scalar.

Parameters

• rel_pos – (npos, 3)-array like with positions in voxel coordinates relative to the voxel centre (betweeen -0.5 and +0.5)

• derivative – computes the derivative in the given dimension (default: no derivative)

Returns

(npos, 2, 2, 2) array which multiplied with the polynomial coefficients gives the scalar value.

gyral_structure.basis.trilinear.grid_to_pos_scalar(rel_pos, derivative=(0, 0, 0))

Given a position relative to the voxel center computes the mapping from the voxel corners to the point.

Parameters

• rel_pos – (npos, 3)-array like with a position in voxel coordinates relative to the voxel centre (betweeen -0.5 and +0.5)

• derivative – computes the derivative in the given dimension (default: no derivative)

Returns

(2, 2, 2) array which multiplied with the scalars in the voxel corers gives the scalar value at the point of interest

gyral_structure.basis.trilinear.grid_to_pos_field(rel_pos, derivative=(0, 0, 0), curl=True)

Computes the mapping from the field potential at the corners to the position of interest.

Parameters

• rel_pos – (npos, 3) array with positions relative to the voxel center of point of interest (between -0.5 and 0.5)

• derivative – computes the derivative in the given dimension (default: no derivative)

• curl – if True computes the fiber field rather than the field potential at the point of interest

Returns

(npos, 3, 2, 2, 2, 3) array which when multiplied with the 3D scalar field in the 2x2x2 neighbourhood gives the field at the point of interest The vector potential at the voxel corners should be represented as a (2, 2, 2, 3) array

gyral_structure.basis.unbounded module

Provides several basis functions:

• ChargeDistribution

• Gradient

• Fourier

• Polynomial

class gyral_structure.basis.unbounded.Fourier(nparams, ndim=3, nfreq=8.0, basefreq=1.0)

Bases: gyral_structure.basis.core.BasisFunc

A set of Fourier basis functions in ndim space.

potential = True

scalar = True

compute = ('matrix', 'numba', 'matrix_cuda')

__init__(nparams, ndim=3, nfreq=8.0, basefreq=1.0)

Creates a new Fourier basis function

Parameters

• nparams – Number of parameters for each of the potential field components

• ndim – Number of physical dimensions

• nfreq – Maximum frequency along any dimension is nfreq times basefreq

• basefreq – lowest frequency

get_partial_mat(positions, derivative=-1)

Computes the matrix used to map parameters to the potential scalar field at the given positions

Parameters

• positions – (npos, 3) array with the positions of interest (shift and scaling have already been applied)

• derivative – which spatial derivative to compute (default: no derivative)

Returns

(npos, nparams // 3) array mapping the parameters to the positions

get_func(positions, method)

Creates a function that evaluates the requested field at the positions.

Parameters

• positions – (npos, ndim) array defining where the result should be computed

• method – which algorithm to use in the computation

Returns

function that maps (input, inverse, derivative) -> output (field derivative to parameter derivative if inverse else param to field)

Raises

NotImplementedError if the basis function is not defined for a specific method

to_hdf5(group: h5py._hl.group.Group)

Stores the radial basis function in the provided HDF5 group

Parameters

group – HDF5 group, where the basis function will be stored

classmethod from_hdf5(group: h5py._hl.group.Group)

Reads the radial basis function from the provided HDF5 group

Parameters

group – HDF5 group, where the basis function is stored

Returns

newly read radial basis function

class gyral_structure.basis.unbounded.ChargeDistribution(charge_pos, charge_size=0.5)

Bases: gyral_structure.basis.core.BasisFunc

Computes the expected fiber orientations given a charge distribution.

Uses the electrostatic model (1/r^2 in 3D and 1/r in 2D) The parameters set how many streamlines terminate at every location

potential = False

scalar = False

compute = ('cuda', 'numba', 'matrix', 'matrix_cuda')

__init__(charge_pos, charge_size=0.5)

Creates a new charge distribution that can be used to allow streamlines to terminate.

Parameters

• charge_pos – Positions where the streamlines terminate (i.e. the locations of the charges)

• charge_size – radius of sphere over which the streamlines uniformly terminate

get_partial_mat(positions)

Computes the matrix used to map parameters to the vector field at the given positions

Parameters

positions – (npos, 3) array with the positions of interest (shift and scaling have already been applied)

Returns

(npos * 3, nparams) array mapping the positions to the parameters

get_func(positions, method)

Creates a function that evaluates the requested field at the positions.

Parameters

• positions – (npos, ndim) array defining where the result should be computed

• method – which algorithm to use in the computation

Returns

function that maps (input, inverse, derivative) -> output (field derivative to parameter derivative if inverse else param to field)

Raises

NotImplementedError if the basis function is not defined for a specific method

to_hdf5(group: h5py._hl.group.Group)

Stores the radial basis function in the provided HDF5 group

Parameters

group – HDF5 group, where the basis function will be stored

classmethod from_hdf5(group: h5py._hl.group.Group)

Reads the radial basis function from the provided HDF5 group

Parameters

group – HDF5 group, where the basis function is stored

Returns

newly read radial basis function

class gyral_structure.basis.unbounded.Polynomial(order, shift=0, scaling=1, ndim=3)

Bases: gyral_structure.basis.core.BasisFunc

Describes the vector potential as a polynomial

potential = True

scalar = True

compute = ('matrix', 'numba', 'matrix_cuda')

__init__(order, shift=0, scaling=1, ndim=3)

Sets up a polynomial basis function with the required polynomial order

Parameters

• order – order of the polynomial

• shift – (ndim, ) array subtracted from positions before calling use_func

• scaling – (ndim, ) array with which the positions are scaled before calling use_func

• ndim – number of spatial dimensions

poly_coeffs()

Computes all the polynomial coefficients.

Returns

(scalar_nparams, ndim) array with the polynomial coefficients

property nparams

Returns the number of parameters needed to describe the field.

get_partial_mat(positions, derivative=-1)

Maps the parameters describing a scalar field to the positions.

Parameters

• positions – (npos, 3) array with the positions of interest (scaling & shift have already been applied to these positions)

• derivative – which spatial derivative to compute (default: no derivative)

Returns

(npos, scalar_nparams) array mapping the parameters to the positions

get_func(positions, method)

Creates a function that evaluates the requested scalar potential at the positions.

Parameters

• positions – (npos, ndim) array defining where the result should be computed

• method – which algorithm to use in the computation

Returns

function that maps (input, inverse, derivative) -> output (derfield to derparam if inverse else param to field)

to_hdf5(group: h5py._hl.group.Group)

Stores the radial basis function in the provided HDF5 group

Parameters

group – HDF5 group, where the basis function will be stored

classmethod from_hdf5(group: h5py._hl.group.Group)

Reads the radial basis function from the provided HDF5 group

Parameters

group – HDF5 group, where the basis function is stored

Returns

newly read radial basis function

class gyral_structure.basis.unbounded.Gradient(order=1, ndim=3, method=())

Bases: gyral_structure.basis.core.BasisFunc

Produces a gradient vector field

The field is given by f = a + b x where a is a vector of length ndim, b is a (ndim x ndim) matrix To get a divergence-free field the trace of b should be zero

potential = False

scalar = False

compute = ('numpy', 'matrix', 'cuda', 'matrix_cuda')

__init__(order=1, ndim=3, method=())

Produces a basic vector field consisting of a constant field + a gradient (latter only if order is 1)

Parameters

• order – if 0 only include a constant field term, otherwise also include a gradient

• ndim – number of dimensions

• method – sorted list of preferred evaluation methods

property nparams

property order

property ndim

int(x=0) -> integer int(x, base=10) -> integer

Convert a number or string to an integer, or return 0 if no arguments are given. If x is a number, return x.__int__(). For floating point numbers, this truncates towards zero.

If x is not a number or if base is given, then x must be a string, bytes, or bytearray instance representing an integer literal in the given base. The literal can be preceded by ‘+’ or ‘-‘ and be surrounded by whitespace. The base defaults to 10. Valid bases are 0 and 2-36. Base 0 means to interpret the base from the string as an integer literal. >>> int(‘0b100’, base=0) 4

get_a(params)

invert_a(dparams)

get_b(params)

invert_b(dparams)

get_partial_mat(positions)

Computes the matrix used to map parameters to the vector field at the given positions

Parameters

positions – (npos, 3) array with the positions of interest

Returns

(npos * 3, nparams) array mapping the parameters to the positions

get_func(positions, method)

Creates a function that evaluates the requested scalar potential at the positions.

Parameters

• positions – (npos, ndim) array defining where the result should be computed

• method – which algorithm to use in the computation

Returns

function that maps (input, inverse, derivative) -> output (derfield to derparam if inverse else param to field)

to_hdf5(group: h5py._hl.group.Group)

Stores the radial basis function in the provided HDF5 group

Parameters

group – HDF5 group, where the basis function will be stored

classmethod from_hdf5(group: h5py._hl.group.Group)

Reads the radial basis function from the provided HDF5 group

Parameters

group – HDF5 group, where the basis function is stored

Returns

newly read radial basis function