mindboggle.features package

Submodules

mindboggle.features.folds module

Functions to extract surface folds.

Authors:

Copyright 2016, Mindboggle team (http://mindboggle.info), Apache v2.0 License

extract_folds(depth_file, depth_threshold=2, min_fold_size=50, save_file=False, output_file='', background_value=-1, verbose=False)

Use depth threshold to extract folds from a triangular surface mesh.

A fold is a group of connected, deep vertices. To extract folds, a depth threshold is used to segment deep vertices of the surface mesh. We have observed in the histograms of travel depth measures of cortical surfaces that there is a rapidly decreasing distribution of low depth values (corresponding to the outer surface, or gyral crowns) with a long tail of higher depth values (corresponding to the folds).

The find_depth_threshold function therefore computes a histogram of travel depth measures, smooths the histogram’s bin values, convolves to compute slopes, and finds the depth value for the first bin with zero slope. The extract_folds function uses this depth value, segments deep vertices, and removes extremely small folds (empirically set at 50 vertices or fewer out of a total mesh size of over 100,000 vertices).

Steps ::
  1. Segment deep vertices as an initial set of folds.
  2. Remove small folds.
  3. Renumber folds.
Note ::
Removed option: Find and fill holes in the folds: Folds could have holes in areas shallower than the depth threshold. Calling fill_holes() could accidentally include very shallow areas (in an annulus-shaped fold, for example). However, we could include the argument exclude_range to check for any values from zero to min_hole_depth; holes would not be filled if they were to contain values within this range.
Parameters:
  • depth_file (string) – surface mesh file in VTK format with faces and depth scalar values
  • depth_threshold (float) – threshold defining the minimum depth for vertices to be in a fold
  • min_fold_size (integer) – minimum fold size (number of vertices)
  • save_file (bool) – save output VTK file?
  • output_file (string) – name of output file in VTK format
  • background_value (integer or float) – background value
  • verbose (bool) – print statements?
Returns:

  • folds (list of integers) – fold numbers for all vertices (-1 for non-fold vertices)
  • n_folds (int) – number of folds
  • folds_file (string (if save_file)) – name of output VTK file with fold IDs (-1 for non-fold vertices)

Examples

>>> from mindboggle.features.folds import extract_folds
>>> from mindboggle.mio.fetch_data import prep_tests
>>> urls, fetch_data = prep_tests()
>>> depth_file = fetch_data(urls['left_travel_depth'], '', '.vtk')
>>> depth_threshold = 2.36089
>>> min_fold_size = 50
>>> save_file = True
>>> output_file = 'extract_folds.vtk'
>>> background_value = -1
>>> verbose = False
>>> folds, n_folds, folds_file = extract_folds(depth_file,
...     depth_threshold, min_fold_size, save_file, output_file,
...     background_value, verbose)
>>> n_folds
33
>>> lens = [len([x for x in folds if x == y]) for y in range(n_folds)]
>>> lens[0:10]
[726, 67241, 2750, 5799, 1151, 6360, 1001, 505, 228, 198]

View folds (skip test):

>>> from mindboggle.mio.plots import plot_surfaces # doctest: +SKIP
>>> plot_surfaces('extract_folds.vtk') # doctest: +SKIP

View folds without background (skip test):

>>> from mindboggle.mio.plots import plot_surfaces # doctest: +SKIP
>>> from mindboggle.mio.vtks import rewrite_scalars # doctest: +SKIP
>>> rewrite_scalars(depth_file, 'extract_folds_no_background.vtk', folds,
...     'just_folds', folds, -1) # doctest: +SKIP
>>> plot_surfaces('extract_folds_no_background.vtk') # doctest: +SKIP
find_depth_threshold(depth_file, min_vertices=10000, verbose=False)

Find depth threshold to extract folds from a triangular surface mesh.

Steps ::
  1. Compute histogram of depth measures.
  2. Define a depth threshold and find the deepest vertices. To extract an initial set of deep vertices from the surface mesh, we anticipate that there will be a rapidly decreasing distribution of low depth values (on the outer surface) with a long tail of higher depth values (in the folds), so we smooth the histogram’s bin values, convolve to compute slopes, and find the depth value for the first bin with slope = 0. This is our threshold.
Parameters:
  • depth_file (string) – surface mesh file in VTK format with faces and depth scalar values
  • min_vertices (integer) – minimum number of vertices
  • verbose (bool) – print statements?
Returns:

  • depth_threshold (float) – threshold defining the minimum depth for vertices to be in a fold
  • bins (list of integers) – histogram bins: each is the number of vertices within a range of depth values
  • bin_edges (list of floats) – histogram bin edge values defining the bin ranges of depth values

Examples

>>> import numpy as np
>>> from mindboggle.features.folds import find_depth_threshold
>>> from mindboggle.mio.fetch_data import prep_tests
>>> urls, fetch_data = prep_tests()
>>> depth_file = fetch_data(urls['left_travel_depth'], '', '.vtk')
>>> min_vertices = 10000
>>> verbose = False
>>> depth_threshold, bins, bin_edges = find_depth_threshold(depth_file,
...     min_vertices, verbose)
>>> np.float("{0:.{1}f}".format(depth_threshold, 5))
2.36089

View threshold histogram plots (skip test):

>>> def vis():
...     import numpy as np
...     import pylab
...     from scipy.ndimage.filters import gaussian_filter1d
...     from mindboggle.mio.vtks import read_scalars
...     # Plot histogram and depth threshold:
...     depths, name = read_scalars(depth_file)
...     nbins = np.round(len(depths) / 100.0)
...     a,b,c = pylab.hist(depths, bins=nbins)
...     pylab.plot(depth_threshold * np.ones((100,1)),
...                np.linspace(0, max(bins), 100), 'r.')
...     pylab.title('Histogram of depth values with threshold')
...     pylab.xlabel('Depth')
...     pylab.ylabel('Number of vertices')
...     pylab.show()
...     # Plot smoothed histogram:
...     bins_smooth = gaussian_filter1d(bins.tolist(), 5)
...     pylab.plot(list(range(len(bins))), bins, '.',
...                list(range(len(bins))), bins_smooth,'-')
...     pylab.title('Smoothed histogram of depth values')
...     pylab.show()
>>> vis() # doctest: +SKIP

mindboggle.features.fundi module

Extract fundus curves from surface mesh patches (folds).

Authors: Arno Klein, 2013-2016 . arno@mindboggle.info . www.binarybottle.com

Copyright 2016, Mindboggle team (http://mindboggle.info), Apache v2.0 License

extract_fundi(folds, curv_file, depth_file, min_separation=10, erode_ratio=0.1, erode_min_size=1, save_file=False, output_file='', background_value=-1, verbose=False)

Extract fundi from folds.

A fundus is a branching curve that runs along the deepest and most highly curved portions of a fold. This function extracts one fundus from each fold by finding the deepest vertices inside the fold, finding endpoints along the edge of the fold, and connecting the former to the latter with tracks that run along deep and curved paths (through vertices with high values of travel depth multiplied by curvature), and a final filtration step.

The deepest vertices are those with values at least two median absolute deviations above the median (non-zero) value, with the higher value chosen if two of the vertices are within (a default of) 10 edges from each other (to reduce the number of possible fundus paths as well as computation time).

To find the endpoints, the find_outer_endpoints function propagates multiple tracks from seed vertices at median depth in the fold through concentric rings toward the fold’s edge, selecting maximal values within each ring, and terminating at candidate endpoints. The final endpoints are those candidates at the end of tracks that have a high median value, with the higher value chosen if two candidate endpoints are within (a default of) 10 edges from each other (otherwise, the resulting fundi can have spurious branching at the fold’s edge).

The connect_points_erosion function connects the deepest fold vertices to the endpoints with a skeleton of 1-vertex-thick curves by erosion. It erodes by iteratively removing simple topological points and endpoints in order of lowest to highest values, where a simple topological point is a vertex that when added to or removed from an object on a surface mesh (such as a fundus curve) does not alter the object’s topology.

Steps ::
  1. Find fundus endpoints (outer anchors) with find_outer_endpoints().
  2. Include inner anchor points.
  3. Connect anchor points using connect_points_erosion(); inner anchors are removed if they result in endpoints.
Note ::
Follow this with segment_by_region() to segment fundi by sulci.
Parameters:
  • folds (numpy array or list of integers) – fold number for each vertex
  • curv_file (string) – surface mesh file in VTK format with mean curvature values
  • depth_file (string) – surface mesh file in VTK format with rescaled depth values
  • likelihoods (list of integers) – fundus likelihood value for each vertex
  • min_separation (integer) – minimum number of edges between inner/outer anchor points
  • erode_ratio (float) – fraction of indices to test for removal at each iteration in connect_points_erosion()
  • save_file (bool) – save output VTK file?
  • output_file (string) – output VTK file
  • background_value (integer or float) – background value
  • verbose (bool) – print statements?
Returns:

  • fundus_per_fold (list of integers) – fundus numbers for all vertices, labeled by fold (-1 for non-fundus vertices)
  • n_fundi_in_folds (integer) – number of fundi
  • fundus_per_fold_file (string (if save_file)) – output VTK file with fundus numbers (-1 for non-fundus vertices)

Examples

>>> # Extract fundus from one or more folds:
>>> import numpy as np
>>> from mindboggle.mio.vtks import read_scalars
>>> from mindboggle.features.fundi import extract_fundi
>>> from mindboggle.mio.fetch_data import prep_tests
>>> urls, fetch_data = prep_tests()
>>> curv_file = fetch_data(urls['left_mean_curvature'], '', '.vtk')
>>> depth_file = fetch_data(urls['left_travel_depth'], '', '.vtk')
>>> folds_file = fetch_data(urls['left_folds'], '', '.vtk')
>>> folds, name = read_scalars(folds_file, True, True)
>>> # Limit number of folds to speed up the test:
>>> limit_folds = True
>>> if limit_folds:
...     fold_numbers = [4] #[4, 6]
...     i0 = [i for i,x in enumerate(folds) if x not in fold_numbers]
...     folds[i0] = -1
>>> min_separation = 10
>>> erode_ratio = 0.10
>>> erode_min_size = 10
>>> save_file = True
>>> output_file = 'extract_fundi_fold4.vtk'
>>> background_value = -1
>>> verbose = False
>>> o1, o2, fundus_per_fold_file = extract_fundi(folds, curv_file,
...     depth_file, min_separation, erode_ratio, erode_min_size,
...     save_file, output_file, background_value, verbose)
>>> lens = [len([x for x in o1 if x == y])
...         for y in np.unique(o1) if y != background_value]
>>> lens[0:10] # [66, 2914, 100, 363, 73, 331, 59, 30, 1, 14] # (if not limit_folds)
[73]

View result without background (skip test):

>>> from mindboggle.mio.plots import plot_surfaces # doctest: +SKIP
>>> from mindboggle.mio.vtks import rewrite_scalars # doctest: +SKIP
>>> rewrite_scalars(fundus_per_fold_file,
...                 'extract_fundi_fold4_no_background.vtk', o1,
...                 'fundus_per_fold', folds) # doctest: +SKIP
>>> plot_surfaces('extract_fundi_fold4_no_background.vtk') # doctest: +SKIP

mindboggle.features.sulci module

Functions to extract sulci from folds.

Authors:

Copyright 2016, Mindboggle team (http://mindboggle.info), Apache v2.0 License

extract_sulci(labels_file, folds_or_file, hemi, min_boundary=1, sulcus_names=[], save_file=False, output_file='', background_value=-1, verbose=False)

Identify sulci from folds in a brain surface according to a labeling protocol that includes a list of label pairs defining each sulcus.

Since folds are defined as deep, connected areas of a surface, and since folds may be connected to each other in ways that differ across brains, there usually does not exist a one-to-one mapping between folds of one brain and those of another. To address the correspondence problem then, we need to find just those portions of the folds that correspond across brains. To accomplish this, Mindboggle segments folds into sulci, which do have a one-to-one correspondence across non-pathological brains. Mindboggle defines a sulcus as a folded portion of cortex whose opposing banks are labeled with one or more sulcus label pairs in the DKT labeling protocol, where each label pair is unique to one sulcus and represents a boundary between two adjacent gyri, and each vertex has one gyrus label.

This function assigns vertices in a fold to a sulcus in one of two cases. In the first case, vertices whose labels are in only one label pair in the fold are assigned to the label pair’s sulcus if they are connected through similarly labeled vertices to the boundary between the two labels. In the second case, the segment_regions function propagates labels from label borders to vertices whose labels are in multiple label pairs in the fold.

Steps for each fold

1. Remove fold if it has fewer than two labels.
2. Remove fold if its labels do not contain a sulcus label pair.
3. Find vertices with labels that are in only one of the fold's
   label boundary pairs. Assign the vertices the sulcus with the label
   pair if they are connected to the label boundary for that pair.
4. If there are remaining vertices, segment into sets of vertices
   connected to label boundaries, and assign a unique ID to each set.
Parameters:
  • labels_file (string) – file name for surface mesh VTK containing labels for all vertices
  • folds_or_file (numpy array, list or string) – fold number for each vertex / name of VTK file containing fold scalars
  • hemi (string) – hemisphere abbreviation in {‘lh’, ‘rh’} for sulcus labels
  • min_boundary (integer) – minimum number of vertices for a sulcus label boundary segment
  • sulcus_names (list of strings) – names of sulci
  • save_file (bool) – save output VTK file?
  • output_file (string) – name of output file in VTK format
  • background_value (integer or float) – background value
  • verbose (bool) – print statements?
Returns:

  • sulci (list of integers) – sulcus numbers for all vertices (-1 for non-sulcus vertices)
  • n_sulci (integers) – number of sulci
  • sulci_file (string) – output VTK file with sulcus numbers (-1 for non-sulcus vertices)

Examples

>>> # Example 1: Extract sulcus from a fold with one sulcus label pair:
>>> import numpy as np
>>> from mindboggle.features.sulci import extract_sulci
>>> from mindboggle.mio.vtks import read_scalars
>>> from mindboggle.mio.fetch_data import prep_tests
>>> urls, fetch_data = prep_tests()
>>> # Load labels, folds, neighbor lists, and sulcus names and label pairs
>>> labels_file = fetch_data(urls['left_freesurfer_labels'], '', '.vtk')
>>> folds_file = fetch_data(urls['left_folds'], '', '.vtk')
>>> folds_or_file, name = read_scalars(folds_file, True, True)
>>> save_file = True
>>> output_file = 'extract_sulci_fold4_1sulcus.vtk'
>>> background_value = -1
>>> # Limit number of folds to speed up the test:
>>> limit_folds = True
>>> if limit_folds:
...     fold_numbers = [4] #[4, 6]
...     i0 = [i for i,x in enumerate(folds_or_file) if x not in fold_numbers]
...     folds_or_file[i0] = background_value
>>> hemi = 'lh'
>>> min_boundary = 10
>>> sulcus_names = []
>>> verbose = False
>>> sulci, n_sulci, sulci_file = extract_sulci(labels_file, folds_or_file,
...     hemi, min_boundary, sulcus_names, save_file, output_file,
...     background_value, verbose)
>>> n_sulci  # 23 # (if not limit_folds)
1
>>> lens = [len([x for x in sulci if x==y])
...         for y in np.unique(sulci) if y != -1]
>>> lens[0:10]  # [6358, 3288, 7612, 5205, 4414, 6251, 3493, 2566, 4436, 739] # (if not limit_folds)
[1151]

View result without background (skip test):

>>> from mindboggle.mio.plots import plot_surfaces # doctest: +SKIP
>>> from mindboggle.mio.vtks import rewrite_scalars # doctest: +SKIP
>>> output = 'extract_sulci_fold4_1sulcus_no_background.vtk'
>>> rewrite_scalars(sulci_file, output, sulci,
...                 'sulci', sulci) # doctest: +SKIP
>>> plot_surfaces(output) # doctest: +SKIP

Example 2: Extract sulcus from a fold with multiple sulcus label pairs:

>>> folds_or_file, name = read_scalars(folds_file, True, True)
>>> output_file = 'extract_sulci_fold7_2sulci.vtk'
>>> # Limit number of folds to speed up the test:
>>> limit_folds = True
>>> if limit_folds:
...     fold_numbers = [7] #[4, 6]
...     i0 = [i for i,x in enumerate(folds_or_file) if x not in fold_numbers]
...     folds_or_file[i0] = background_value
>>> sulci, n_sulci, sulci_file = extract_sulci(labels_file, folds_or_file,
...     hemi, min_boundary, sulcus_names, save_file, output_file,
...     background_value, verbose)
>>> n_sulci  # 23 # (if not limit_folds)
2
>>> lens = [len([x for x in sulci if x==y])
...         for y in np.unique(sulci) if y != -1]
>>> lens[0:10]  # [6358, 3288, 7612, 5205, 4414, 6251, 3493, 2566, 4436, 739] # (if not limit_folds)
[369, 93]

View result without background (skip test):

>>> from mindboggle.mio.plots import plot_surfaces # doctest: +SKIP
>>> from mindboggle.mio.vtks import rewrite_scalars # doctest: +SKIP
>>> output = 'extract_sulci_fold7_2sulci_no_background.vtk'
>>> rewrite_scalars(sulci_file, output, sulci,
...                 'sulci', sulci) # doctest: +SKIP
>>> plot_surfaces(output) # doctest: +SKIP

Module contents