|
| 1 | +import os |
| 2 | +from argparse import ArgumentParser |
| 3 | + |
| 4 | +import nipype.pipeline.engine as pe |
| 5 | +from niworkflows.interfaces.nibabel import MapLabels |
| 6 | + |
| 7 | +from nibabies.interfaces.mcribs import MCRIBReconAll |
| 8 | + |
| 9 | + |
| 10 | +def _parser(): |
| 11 | + parser = ArgumentParser(description="Test script for MCRIBS surfaces") |
| 12 | + parser.add_argument('subject', help='Subject ID') |
| 13 | + parser.add_argument('t2w', type=os.path.abspath, help='Input T2w (radioisotropic)') |
| 14 | + parser.add_argument( |
| 15 | + 'segmentation', type=os.path.abspath, help='Input anatomical segmentation in T2w space' |
| 16 | + ) |
| 17 | + parser.add_argument( |
| 18 | + '--outdir', type=os.path.abspath, help='Output directory to persist MCRIBS output' |
| 19 | + ) |
| 20 | + parser.add_argument('--nthreads', type=int, help='Number of threads to parallelize tasks') |
| 21 | + return parser |
| 22 | + |
| 23 | + |
| 24 | +def main(argv: list = None): |
| 25 | + pargs = _parser().parse_args(argv) |
| 26 | + |
| 27 | + t2w_file = _check_file(pargs.t2w) |
| 28 | + seg_file = _check_file(pargs.segmentation) |
| 29 | + |
| 30 | + aseg2mcrib = { |
| 31 | + 2: 51, |
| 32 | + 3: 21, |
| 33 | + 4: 49, |
| 34 | + 5: 0, |
| 35 | + 7: 17, |
| 36 | + 8: 17, |
| 37 | + 10: 43, |
| 38 | + 11: 41, |
| 39 | + 12: 47, |
| 40 | + 13: 47, |
| 41 | + 14: 0, |
| 42 | + 15: 0, |
| 43 | + 16: 19, |
| 44 | + 17: 1, |
| 45 | + 18: 3, |
| 46 | + 26: 41, |
| 47 | + 28: 45, |
| 48 | + 31: 49, |
| 49 | + 41: 52, |
| 50 | + 42: 20, |
| 51 | + 43: 50, |
| 52 | + 44: 0, |
| 53 | + 46: 18, |
| 54 | + 47: 18, |
| 55 | + 49: 42, |
| 56 | + 50: 40, |
| 57 | + 51: 46, |
| 58 | + 52: 46, |
| 59 | + 53: 2, |
| 60 | + 54: 4, |
| 61 | + 58: 40, |
| 62 | + 60: 44, |
| 63 | + 63: 50, |
| 64 | + 253: 48, |
| 65 | + } |
| 66 | + map_labels = pe.Node(MapLabels(in_file=seg_file, mappings=aseg2mcrib), name='map_labels') |
| 67 | + |
| 68 | + recon = pe.Node( |
| 69 | + MCRIBReconAll(subject_id=pargs.subject, t2w_file=t2w_file), name='mcribs_recon' |
| 70 | + ) |
| 71 | + if pargs.outdir: |
| 72 | + recon.inputs.outdir = pargs.outdir |
| 73 | + if pargs.nthreads: |
| 74 | + recon.inputs.nthreads = pargs.nthreads |
| 75 | + |
| 76 | + wf = pe.Workflow(f'MRA_{pargs.subject}') |
| 77 | + wf.connect(map_labels, 'out_file', recon, 'segmentation_file') |
| 78 | + wf.run() |
| 79 | + |
| 80 | + |
| 81 | +def _check_file(fl: str) -> str: |
| 82 | + import nibabel as nb |
| 83 | + import numpy as np |
| 84 | + |
| 85 | + img = nb.load(fl) |
| 86 | + if len(img.shape) != 3: |
| 87 | + raise ValueError('Image {fl} is not 3 dimensional.') |
| 88 | + |
| 89 | + voxdims = img.header['pixdim'][1:4] |
| 90 | + if not np.allclose(voxdims, voxdims[1]): |
| 91 | + raise ValueError(f'Image {fl} is not isotropic: {voxdims}.') |
| 92 | + |
| 93 | + ornt = nb.io_orientation(img.affine) |
| 94 | + axcodes = nb.orientations.ornt2axcodes(ornt) |
| 95 | + if ''.join(axcodes) != 'LAS': |
| 96 | + las = nb.orientations.axcodes2ornt('LAS') |
| 97 | + transform = nb.orientations.ornt_transform(ornt, las) |
| 98 | + reornt = img.as_reoriented(transform) |
| 99 | + outfl = os.path.abspath(f'LASornt_{os.path.basename(fl)}') |
| 100 | + print(f'Creating reorientated image {outfl}') |
| 101 | + reornt.to_filename(outfl) |
| 102 | + return outfl |
| 103 | + return fl |
| 104 | + |
| 105 | + |
| 106 | +if __name__ == '__main__': |
| 107 | + main() |
0 commit comments