Tutorial

Installation

Gyral structure can be installed using pip:

pip install git+https://git.fmrib.ox.ac.uk/ndcn0236/mcutils.git
pip install git+https://git.fmrib.ox.ac.uk/ndcn0236/gyral_structure.git

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

Command line interface

After installation several scripts will be made available to model a continuous vector field in the gyral white matter. These scripts will require reconstruction of the surfaces (e.g., from freesurfer) as well as estimates of the main fibre orientation in the gyral white matter (e.g., from DTI). For a preprocessed dataset from the Human Connectome Project these files would for example be available as:

subject=100307
hcpdir=<directory with Human Connectome Project data>
WHITE=${hcpdir}/${subject}/T1w/fsaverage_32k/${subject}.L.white.32k_fs_LR.surf.gii
MID=${hcpdir}/${subject}/T1w/fsaverage_32k/${subject}.L.midthickness.32k_fs_LR.surf.gii  # mid-thickness surface is optional
PIAL=${hcpdir}/${subject}/T1w/fsaverage_32k/${subject}.L.pial.32k_fs_LR.surf.gii
SURF_MASK=${hcpdir}/${subject}/MNINonLinear/fsaverage_32k/${subject}.L.thickness.32k_fs_LR.shape.gii
VOL_REF=${hcpdir}/${subject}/T1w/Diffusion.bedpostX/nodif_brain_mask.nii.gz
DYADS=${hcpdir}/${subject}/T1w/Diffusion.bedpostX/dyads1.nii.gz

We can then model the fibre configuration in the gyral white matter and generate a new surface at the interface between the gyral and deep white matter using:

gs_mask $WHITE $SURF_MASK $VOL_REF gyral_wm_mask.nii.gz
gs_fit gyral_wm_mask.nii.gz $WHITE $PIAL $SURF_MASK $DYADS best_fit.pck.gz -m $MID -ns 1 -nv 1
gs_deform_surface best_fit.pck.gz $WHITE gyral_wm_mask.nii.gz deep_gyral_interface.surf.gii

The first script gs_mask creates a mask of the gyral white matter. The two main variables controlling this script are the width (-w or –width) which sets the maximum distance between the gyral walls between which the white matter is still considered to be gyral and the distance (-d or –dist) which sets the maximum distance from the gyral walls that will be included in the gyral white matter even if there is no opposite gyral wall (i.e., include WM just below the sulcal fundi). This script can be made faster by reducing the number of orientations that will be used to estimate the gyral width (-N), although this will make the width estimates less accurate.

The next script gs_fit does the actual fitting of the vector field to the gyral white matter. This script is described in much more detail below. The field is initialised as the sum of the field from negative charges uniformly distributed on the pial surface with an equal-sized possitive charge in the deep white matter. To this field radial basis functions are added, whose strength and orientation is determined by minimising a cost function. This cost function has four terms:

  • Ensure smooth fibre density across both the white/gray matter boundary and the mid-cortical surface

  • Radial orientation at both the white/gray matter boundary and the mid-cortical surface

  • Alignment with the dyads in the gyral white matter

  • L2 norm on the fibre density in the gyral white matter

The final script gs_deform_surface propagates the white/gray matter boundary through the best-fit vector field to generate a new surface at the interface between the deep and gyral white matter.

Other scripts can be used to evaluate the fibre orientation and density of the best-fit vector field on a volumetric (gs_dyad_vol) or surface (gs_dyad_surf) mask or to do streamline tractography (gs_track).

Deconstructing gs_fit

To become more familiar with the gyral_structure python API we will here take us step-by-step through (a somewhat simplified version of) the gs_fit code.

After defining the command line interface we start by importing the libraries:

import nibabel as nib
import numpy as np
from gyral_structure import basis, cost, request
from mcutils.surface import CorticalMesh, grid

Requests: where to evaluate the cost function

After imports the first step is to define where the vector field will be evaluated:

white_request = request.read_surface(args.white, args.surf_mask, nsplit=args.nsplit_surf)
if args.mid_thickness is not None:
    mid_request = request.read_surface(args.mid_thickness, args.surf_mask, nsplit=args.nsplit_surf)
vol_request = request.read_volume(args.gyral_wm_mask, nsplit=args.nsplit_vol)

The variables white_request and optionally mid_request define where on the surface we want to evaluate the vector field (i.e., an instance of request.VertexRequest). The vector field is evaluated for each triangle on the surface mesh, whose vertices are in the surface mask. The nsplit variable determines at how many points on the triangle the vector field will be evaluated when attempting to estimate the mean vector field. Each triangle will be split into \(({\rm nsplit})^2\) equal-sized sub-triangles with the field for each being evaluated at the sub-triangle’s central point.

Similarly vol_request will be a request.VoxelRequest instance defining for which voxels the field will be evaluated. Each voxel’s field will be evaluated as the average of field at the center of \(({\rm nsplit})^3\) sub-voxels. When setting nsplit it should be kept in mind that if the matrix mapping the parameters to the vector field is not precomputed, a higher nsplit will cause a much longer runtime (if the matrix mapping is precomputed, it will only slow down the preprocessing, which might not significantly alter the total runtime).

Cost: defining the cost function

Once we have defined where the field will be evaluated, we can define the individual terms to the cost function. The simplest of these is:

L2 = cost.dens.L2Density(vol_request)

L2Density is like all cost function terms a sub-class of CostFunc. Each cost function expect as their first term a request defining to which vectors the cost function should be applied (the only exception to this are the cost functions in cost.param, which are evaluated directly on the model parameters rather than the vector field.

Additional terms to further define the cost function might have to be provided as well. For example when defining the term minimizing the offset with the observed dyads, we have to provide those observed dyads as well:

dyads_arr = nib.load(args.dyads).get_fdata()[nib.load(args.gyral_wm_mask).get_fdata() != 0]
affine = nib.load(args.dyads).affine[:3, :3]
if np.linalg.det(affine) > 0:
    # bvecs assume negative determinant, so need to flip around x-axis
    dyads_arr[..., 0] *= -1
mm_dyads_arr = affine.dot(dyads_arr.T).T
mm_dyads_arr /= np.sqrt((mm_dyads_arr ** 2).sum(-1))[:, None]
aligned = cost.orient.Watson(vol_request, field=mm_dyads_arr)

Note that we fit the field in mm rather than scaled voxel space, so we need to convert the dyads from radiological voxel space to mm space.

Before defining the cost function on the surface, we need to load the target density distribution. For each triangular face in the surface, the target number of streamlines is proportional to the size of the wedge covered by that triangle:

target_density = Cortex([white_surf, CorticalMesh.read(args.pial)[surf_mask]]).wedge_volume()

This allows us to define the cost function at the surface:

white_normal = white_surf.normal().T
white_radial = cost.orient.VonMises(white_request, field=white_normal) / white_request.npos
white_density = cost.dens.TargetInner(white_request, field=white_normal,
                                      target_inner=target_density / white_surf.size_faces()) / white_request.npos
if args.mid_thickness is not None:
    mid_surf = CorticalMesh.read(args.mid_thickness)[surf_mask]
    mid_normal = mid_surf.normal().T
    mid_radial = cost.orient.VonMises(mid_request, field=mid_normal) / mid_request.npos
    mid_uniform = cost.dens.TargetInner(mid_request, field=mid_normal,
                                        target_inner=target_density / mid_surf.size_faces()) / mid_request.npos
    radial = 0.5 * (white_radial + mid_radial)
    density = 0.5 * (white_density + mid_uniform)
else:
    radial = white_radial
    density = white_density

If no mid-thickness surface was provided (i.e., args.mid_thickness is None), than the resulting cost functions are simply instances of cost.orient.VonMises (to encourage radiality) and cost.dens.TargetInner (to encourage a mean streamline termination density of 1 per mm$^3$).

If a mid-thickness is provided we define the same cost functions for the request of the mid-thickness surface. These are then combined into a single cost-function by adding them together and multiplying by 0.5 to keep the same total weight, as can be seen in the lines:

radial = 0.5 * (white_radial + mid_radial)
uniform = 0.5 * (white_uniform + mid_uniform)

Once all the individual terms are defined we can add them together with the user-provided weights:

no_dyads = args.wdensity * density + args.wradial * radial + args.wL2 * L2
with_dyads = no_dyads + args.wdyad * aligned

Note that this adding of cost-functions is done before any evaluation of the vector field. It merely provides instructions for the optimisation later on how these cost functions should be combined.

Basis: mapping parameters to the vector field

Now that we have a cost function, we still have to define our basis functions:

white = CorticalMesh.read(args.white)
ref_img = nib.load(args.gyral_wm_mask)
white_dist = grid.signed_distance(white, ref_img.shape, affine=ref_img.affine)
idx_max = np.argmin(white_dist)
voxel_max = np.array(np.unravel_index(idx_max, ref_img.shape)).flatten()
mm_max = ref_img.affine[:3, :3].dot(voxel_max) + ref_img.affine[:3, -1]

masked_pial = CorticalMesh.read(args.pial)[surf_mask]
points_pial = masked_pial.vertices[:, masked_pial.faces].mean(1).T

charge_pial = -target_density

charge_max = -np.atleast_1d(charge_pial.sum())

charge_basis = basis.SumBase([
    basis.ChargeDistribution(points_pial, 0.1),
    basis.ChargeDistribution(mm_max[None, :], 0.1)
])

This code block defines the initial vector field estimate. As our cost function expects a target_density streamlines crossing the surface, we can create a rough approximation of this by distributing uniform charges of -target_density along the pial surface. This basis function is defined by the line basis.ChargeDistribution(points_pial, 0.1) which defines where the charges are (i.e., points_pial) and the spatial extent of the charges (i.e., 0.1 mm).

Our final basis function charge_basis combines this basis function with one defined by a single charge at the centre of the brain (i.e., the point furthest away from the white/gray matter boundary). For both of these parts of the basis function we fix the parameters to their pre-computed charges:

charge_basis.fix(0, charge_pial)
charge_basis.fix(1, charge_max)

This ensures that the charges are kept fixed and not altered during the optimisation (which would be very slow).

Finally we define the radial basis functions, which is the part of the basis function we will actually want to optimise during the fitting:

def get_rbf(resolution=3., spacing=1./3.):
    """
    Defines the radial basis functions

    :param resolution: radius of the basis functions
    :param spacing: spacing between basis functions relative to the resolution (best between 1/2 and 1/3)
    :return: radial basis functions
    """
    req = request.MultRequest([white_request, vol_request])
    bb = req.bounding_box()
    centroids = basis.hcp_packing(bb, distance=resolution * spacing)
    use_centroids = req.min_dist(centroids) < resolution / 2.
    return basis.RadialBasis(basis.wendland(smoothness=2), centroids[use_centroids], resolution)

The first two lines are simply meant to get a bounding-box around the request, so that we can spread the centroids of the radial basis functions throughout the region, where we actually want to evaluate the vector field. These centroids are computed using basis.hcp_packing. From the centroids throughout this bounding box we only include those that are within \(0.5 * size_{RBF}\) of the points where we evaluate the cost function (which can be usefully evaluated using FieldRequest.min_dist).

Finally, we return basis.RadialBasis instance that maps the parameters to the vector field. The shape of the radial function can be set using equations from basis.wendland or basis.buhmann.

Optimising basis function parameters to minimise cost function

Because the solution space of one of the terms in the cost function (i.e., alignment with the dyads) is bimodal, we will first do a rough fit without this term in the cost function:

init_basis = basis.SumBase([charge_basis, get_rbf(args.first_rbf)])
fitter = no_dyads.get_evaluator(init_basis)
res = fitter.fit(np.zeros(init_basis.nparams))
print('Initial fitting:')
print(res)
init_basis.fix(1, res.x)

We combine our initial field configuration defined by charge_basis with radial basis functions with radius of 20 mm to do this initial fit. This new basis function init_basis is then combined with the cost function to create fitter. fitter is an instance of the CostEvaluator class, which contains methods to actually fit the vector field.

In this case we simply use the fit method, which can use any optimiser available for scipy.optimize.minimize() to fit the basis function parameters to minimise the cost function. Both the derivative and the hessian multiplied with an arbitrary vector (i.e., hessp) can be efficiently evaluated for all cost functions, so we highly recommend to use an optimisation strategy that uses this additional information (note that typically the number of parameters is too large to evaluate the full hessian).

During the next optimisation stage, we will add a new set of radial basis functions with a smaller radius, so we choose to fix the parameters in this set in the final line (init_basis.fix(1, res.x)). In principle, it is not strictly necessary to do so (we can optimise the parameters for multiple parts of the basis function at the same time).

final_basis = basis.SumBase([init_basis, get_rbf(args.second_rbf)])
fitter = with_dyads.get_evaluator(final_basis)
res = fitter.fit(np.zeros(final_basis.nparams))
print('Final fitting:')
print(res)
final_basis.fix(1, res.x)
final_basis.save(args.out)

The final fitting block is very similar to the previous one, except now we optimise radial basis functions with a much smaller radius (7 instead of 20 mm) and we employ the full cost function with the term aligning the vector field with the dyads (i.e., with_dyads instead of no_dyads).

The final basis function is saved to the output HDF5 file chosen by the user (args.out), so that we can evaluate it during the tractography in the later scripts. Note that again we could have chosen not to fix the basis function with our best-fit parameters, in which case we would have had to pass those parameters as a second term in the BasisFunc.save method.