Skip to content

Commit

Permalink
FIX: Revert to FAST for tissue probability maps
Browse files Browse the repository at this point in the history
  • Loading branch information
effigies committed Aug 30, 2021
1 parent d160455 commit 4e1ea40
Showing 1 changed file with 24 additions and 93 deletions.
117 changes: 24 additions & 93 deletions smriprep/workflows/anatomical.py
Original file line number Diff line number Diff line change
Expand Up @@ -475,38 +475,39 @@ def _check_img(img):
])
# fmt:on

if not freesurfer: # Flag --fs-no-reconall is set - return
# Brain tissue segmentation - FAST produces: 0 (bg), 1 (wm), 2 (csf), 3 (gm)
t1w_dseg = pe.Node(
fsl.FAST(segments=True, no_bias=True, probability_maps=True),
name="t1w_dseg",
mem_gb=3,
)
lut_t1w_dseg.inputs.lut = (0, 3, 1, 2) # Maps: 0 -> 0, 3 -> 1, 1 -> 2, 2 -> 3.
fast2bids = pe.Node(
niu.Function(function=_probseg_fast2bids),
name="fast2bids",
run_without_submitting=True,
)
# XXX Keeping FAST separate so that it's easier to swap in ANTs or FreeSurfer

# Brain tissue segmentation - FAST produces: 0 (bg), 1 (wm), 2 (csf), 3 (gm)
t1w_dseg = pe.Node(
fsl.FAST(segments=True, no_bias=True, probability_maps=True),
name="t1w_dseg",
mem_gb=3,
)
lut_t1w_dseg.inputs.lut = (0, 3, 1, 2) # Maps: 0 -> 0, 3 -> 1, 1 -> 2, 2 -> 3.
fast2bids = pe.Node(
niu.Function(function=_probseg_fast2bids),
name="fast2bids",
run_without_submitting=True,
)

# fmt:off
workflow.connect([
(buffernode, t1w_dseg, [('t1w_brain', 'in_files')]),
(t1w_dseg, lut_t1w_dseg, [('partial_volume_map', 'in_dseg')]),
(t1w_dseg, fast2bids, [('partial_volume_files', 'inlist')]),
(fast2bids, anat_norm_wf, [('out', 'inputnode.moving_tpms')]),
(fast2bids, outputnode, [('out', 't1w_tpms')]),
])
# fmt:on
if not freesurfer: # Flag --fs-no-reconall is set - return
# fmt:off
workflow.connect([
(brain_extraction_wf, buffernode, [
(('outputnode.out_file', _pop), 't1w_brain'),
('outputnode.out_mask', 't1w_mask')]),
(buffernode, t1w_dseg, [('t1w_brain', 'in_files')]),
(t1w_dseg, lut_t1w_dseg, [('partial_volume_map', 'in_dseg')]),
(t1w_dseg, fast2bids, [('partial_volume_files', 'inlist')]),
(fast2bids, anat_norm_wf, [('out', 'inputnode.moving_tpms')]),
(fast2bids, outputnode, [('out', 't1w_tpms')]),
])
# fmt:on
return workflow

# Map FS' aseg labels onto three-tissue segmentation
lut_t1w_dseg.inputs.lut = _aseg_to_three()
split_seg = pe.Node(niu.Function(function=_split_segments), name="split_seg")

# check for older IsRunning files and remove accordingly
fs_isrunning = pe.Node(
niu.Function(function=_fs_isRunning), overwrite=True, name="fs_isrunning"
Expand Down Expand Up @@ -537,12 +538,6 @@ def _check_img(img):
(('outputnode.bias_corrected', _pop), 'in_file')]),
(surface_recon_wf, applyrefined, [
('outputnode.out_brainmask', 'mask_file')]),
(surface_recon_wf, lut_t1w_dseg, [
('outputnode.out_aseg', 'in_dseg')]),
(lut_t1w_dseg, split_seg, [('out', 'in_file')]),
(split_seg, anat_norm_wf, [
('out', 'inputnode.moving_tpms')]),
(split_seg, outputnode, [('out', 't1w_tpms')]),
(surface_recon_wf, outputnode, [
('outputnode.subjects_dir', 'subjects_dir'),
('outputnode.subject_id', 'subject_id'),
Expand Down Expand Up @@ -742,70 +737,6 @@ def _pop(inlist):
return inlist


def _aseg_to_three():
"""
Map FreeSurfer's segmentation onto a brain (3-)tissue segmentation.
This function generates an index of 255+0 labels and maps them into zero (bg),
1 (GM), 2 (WM), or 3 (CSF). The new values are set according to BIDS-Derivatives.
Then the index is populated (e.g., label 3 in the original segmentation maps to label
1 in the output).
The `aseg lookup table
<https://github.com/freesurfer/freesurfer/blob/2beb96c6099d96508246c14a24136863124566a3/distribution/ASegStatsLUT.txt>`__
is available in the FreeSurfer source.
"""
import numpy as np

# Base struct
aseg_lut = np.zeros((256,), dtype="int")
# GM
aseg_lut[3] = 1
aseg_lut[8:14] = 1
aseg_lut[17:21] = 1
aseg_lut[26:40] = 1
aseg_lut[42] = 1
aseg_lut[47:73] = 1

# CSF
aseg_lut[4:6] = 3
aseg_lut[14:16] = 3
aseg_lut[24] = 3
aseg_lut[43:45] = 3
aseg_lut[72] = 3

# WM
aseg_lut[2] = 2
aseg_lut[7] = 2
aseg_lut[16] = 2
aseg_lut[28] = 2
aseg_lut[41] = 2
aseg_lut[46] = 2
aseg_lut[60] = 2
aseg_lut[77:80] = 2
aseg_lut[250:256] = 2
return tuple(aseg_lut)


def _split_segments(in_file):
from pathlib import Path
import numpy as np
import nibabel as nb

segimg = nb.load(in_file)
data = np.int16(segimg.dataobj)
hdr = segimg.header.copy()
hdr.set_data_dtype("uint8")

out_files = []
for i, label in enumerate(("GM", "WM", "CSF"), 1):
out_fname = str(Path.cwd() / f"aseg_label-{label}_mask.nii.gz")
segimg.__class__(data == i, segimg.affine, hdr).to_filename(out_fname)
out_files.append(out_fname)

return out_files


def _probseg_fast2bids(inlist):
"""Reorder a list of probseg maps from FAST (CSF, WM, GM) to BIDS (GM, WM, CSF)."""
return (inlist[1], inlist[2], inlist[0])

0 comments on commit 4e1ea40

Please sign in to comment.