Commit 3a64ff2a authored by Michiel Cottaar's avatar Michiel Cottaar
Browse files

RF: split gs_dyad in gs_dyad_surf and gs_dyad_vol

parent 608bfc11
Pipeline #5228 failed with stage
in 8 minutes and 11 seconds
......@@ -27,11 +27,18 @@ gs_deform_surface: Moves white/gray matter boundary to deep white matter
.. literalinclude:: ../gyral_structure/scripts/gs_deform_surface
:language: python
.. _gs_dyad:
.. _gs_dyad_vol:
gs_dyad: Evaluates vector field fibre density and orientation
-------------------------------------------------------------
.. literalinclude:: ../gyral_structure/scripts/gs_dyad
gs_dyad_vol: Evaluates vector field fibre density and orientation in a volume
-----------------------------------------------------------------------------
.. literalinclude:: ../gyral_structure/scripts/gs_dyad_vol
:language: python
.. _gs_dyad_surf:
gs_dyad_surf: Evaluates vector field fibre density and orientation on a surface
-------------------------------------------------------------------------------
.. literalinclude:: ../gyral_structure/scripts/gs_dyad_surf
:language: python
.. _gs_track:
......
......@@ -63,8 +63,8 @@ The final script :ref:`gs_deform_surface <gs_deform_surface>` propagates the whi
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 or surface mask (:ref:`gs_dyad <gs_dyad>`) or to do streamline tractography from
arbitrary seed points (:ref:`gs_track <gs_track>`).
on a volumetric (:ref:`gs_dyad_vol <gs_dyad_vol>`) or surface (:ref:`gs_dyad_surf <gs_dyad_surf>`) mask
or to do streamline tractography (:ref:`gs_track <gs_track>`).
.. _dec-gs_fit:
......
#!/usr/bin/env python
import argparse
parser = argparse.ArgumentParser('Computes the fibre orientation & density in a volume image')
parser.add_argument('field', help='Input HDF5 file with the basis function and best-fit array')
parser.add_argument('mask', help='reference volume/surface on which the dyads will be computed for all non-zero voxels/vertices')
parser.add_argument('dyad', help='Output volume/surface with fibre orientations')
parser.add_argument('-d', '--density', help='Output volume with fibre densities')
parser.add_argument('--nsplit', type=int, default=1, help='number of sub-voxels/sub-triangles to split each element into along each dimension (default: 1)')
parser.add_argument('--no-flipx', action='store_true', help="don't flip the dyads in the x-direction")
args = parser.parse_args()
import nibabel
from gyral_structure import basis, request, utils
from mcutils.surface import CorticalMesh
import scipy as sp
import cifti
basis_func, best_fit_arr = basis.read(args.field)
img = nibabel.load(args.mask)
if isinstance(img, nibabel.Nifti1Image):
mask = img.get_data() != 0
req = request.read_volume(img, nsplit=args.nsplit)
else:
req = request.read_surface(img, nsplit=args.nsplit)
mask = sp.ones(req.npos, dtype='bool')
field = req.get_evaluator(basis_func)(best_fit_arr)
density = sp.zeros(mask.shape)
density[mask] = sp.sqrt((field ** 2).sum(-1))
dyad = sp.zeros(mask.shape + (3, ))
if isinstance(img, nibabel.Nifti1Image):
dyad[mask] = utils.dyad_from_mm(field, img.affine)
else:
dyad[mask] = field
if isinstance(img, nibabel.Nifti1Image):
nibabel.Nifti1Image(dyad, img.affine).to_filename(args.dyad)
if args.density is not None:
nibabel.Nifti1Image(density, img.affine).to_filename(args.density)
else:
surf = CorticalMesh.read(img)
normals = surf.normal()
crossing = (normals.T * field).sum(-1)
arr_triangle = sp.stack((crossing, density, crossing / density, dyad[:, 0], dyad[:, 1], dyad[:, 2]), axis=0)
arr_vertices = sp.zeros((arr_triangle.shape[0], surf.nvertices), dtype='f4')
for idx in range(arr_triangle.shape[0]):
conn, val = sp.broadcast_arrays(surf.faces, arr_triangle[idx][None, :])
nconn = sp.bincount(conn.flatten(), minlength=surf.nvertices)
arr_vertices[idx] = sp.bincount(conn.flatten(), val.flatten(), minlength=surf.nvertices) / nconn
bm = cifti.BrainModel.from_mask(sp.ones(surf.nvertices, dtype='bool'),
name=img.darrays[0].metadata['AnatomicalStructurePrimary'])
cifti.write(args.dyad, arr_vertices, (cifti.Scalar.from_names(['density (crossing)', 'density (total)', 'radiality', 'x', 'y', 'z']), bm))
#!/usr/bin/env python
import argparse
parser = argparse.ArgumentParser('Computes the fibre orientation & density on a surface')
parser.add_argument('field', help='Input HDF5 file with the basis function and best-fit array')
parser.add_argument('surface', help='reference surface on which the dyads will be computed')
parser.add_argument('output', help='Output CIFTI (.dscalar.nii) file with densities and dyads')
parser.add_argument('--nsplit', type=int, default=1, help='number of sub-voxels/sub-triangles to split each element into along each dimension (default: 1)')
parser.add_argument('-m', '--mask', help='GIFTI filename with surface mask (all non-zero vertices will be included)')
args = parser.parse_args()
import nibabel
from gyral_structure import basis, request, utils
from mcutils.surface import CorticalMesh
import scipy as sp
import cifti
basis_func, best_fit_arr = basis.read(args.field)
img = nibabel.load(args.surface)
req = request.read_surface(img, args.mask, nsplit=args.nsplit)
field = req.get_evaluator(basis_func)(best_fit_arr)
surf = CorticalMesh.read(img)
density = sp.zeros(field.shape[0])
density = sp.sqrt((field ** 2).sum(-1))
normals = surf.normal()
crossing = (normals.T * field).sum(-1)
arr_triangle = sp.stack((crossing, density, crossing / density, field[:, 0], field[:, 1], field[:, 2]), axis=0)
arr_vertices = sp.zeros((arr_triangle.shape[0], surf.nvertices), dtype='f4')
for idx in range(arr_triangle.shape[0]):
conn, val = sp.broadcast_arrays(surf.faces, arr_triangle[idx][None, :])
nconn = sp.bincount(conn.flatten(), minlength=surf.nvertices)
arr_vertices[idx] = sp.bincount(conn.flatten(), val.flatten(), minlength=surf.nvertices) / nconn
bm = cifti.BrainModel.from_mask(sp.ones(surf.nvertices, dtype='bool'),
name=img.darrays[0].metadata['AnatomicalStructurePrimary'])
cifti.write(args.dyad, arr_vertices, (cifti.Scalar.from_names(['density (crossing)', 'density (total)', 'radiality', 'x', 'y', 'z']), bm))
#!/usr/bin/env python
import argparse
parser = argparse.ArgumentParser('Computes the fibre orientation and optionally density in a volume image')
parser.add_argument('field', help='Input HDF5 file with the basis function and best-fit array')
parser.add_argument('mask', help='reference volume on which the dyads will be computed for all non-zero voxels')
parser.add_argument('dyad', help='Output volume filename with fibre orientations')
parser.add_argument('-d', '--density', help='Output volume filename with fibre densities')
parser.add_argument('--nsplit', type=int, default=1, help='number of sub-voxels/sub-triangles to split each element into along each dimension (default: 1)')
args = parser.parse_args()
import nibabel
from gyral_structure import basis, request, utils
import scipy as sp
basis_func, best_fit_arr = basis.read(args.field)
img = nibabel.load(args.mask)
mask = img.get_data() != 0
req = request.read_volume(img, nsplit=args.nsplit)
field = req.get_evaluator(basis_func)(best_fit_arr)
density = sp.zeros(mask.shape)
density[mask] = sp.sqrt((field ** 2).sum(-1))
dyad = sp.zeros(mask.shape + (3, ))
dyad[mask] = utils.dyad_from_mm(field, img.affine)
nibabel.Nifti1Image(dyad, img.affine).to_filename(args.dyad)
if args.density is not None:
nibabel.Nifti1Image(density, img.affine).to_filename(args.density)
......@@ -46,7 +46,8 @@ setup(name='gyral_structure',
install_requires=['nibabel', 'scipy', 'mcutils', 'sympy', 'h5py', 'numexpr'],
scripts=[
'gyral_structure/scripts/gs_track',
'gyral_structure/scripts/gs_dyad',
'gyral_structure/scripts/gs_dyad_vol',
'gyral_structure/scripts/gs_dyad_surf',
'gyral_structure/scripts/gs_deform_surface',
'gyral_structure/scripts/gs_mask',
'gyral_structure/scripts/gs_fit',
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment