From ace3e4efd2266f60f8abe4dd125d2c1edda939aa Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 20 Aug 2020 16:21:59 -0400 Subject: [PATCH 1/7] ENH: Add nib-roi command to crop (maybe flip) axes --- nibabel/cmdline/roi.py | 64 ++++++++++++++++++++++++++++++++++++++++++ setup.cfg | 1 + 2 files changed, 65 insertions(+) create mode 100644 nibabel/cmdline/roi.py diff --git a/nibabel/cmdline/roi.py b/nibabel/cmdline/roi.py new file mode 100644 index 0000000000..485548c08e --- /dev/null +++ b/nibabel/cmdline/roi.py @@ -0,0 +1,64 @@ +import argparse +import nibabel as nb + + +def lossless_slice(img, slicers): + if not nb.imageclasses.spatial_axes_first(img): + raise ValueError("Cannot slice an image that is not known to have spatial axes first") + + roi_img = img.__class__( + img.dataobj._get_unscaled(slicers), + affine=img.slicer.slice_affine(slicers), + header=img.header) + roi_img.header.set_slope_inter(img.dataobj.slope, img.dataobj.inter) + return roi_img + + +def parse_slice(crop, allow_step=True): + if crop is None: + return slice(None) + start, stop, *extra = [int(val) if val else None for val in crop.split(":")] + if len(extra) > 1: + raise ValueError(f"Cannot parse specification: {crop}") + if extra and not allow_step: + raise ValueError(f"Step entry not permitted: {crop}") + + step = extra[0] if extra else None + if step not in (1, -1, None): + raise ValueError(f"Downsampling is not supported: {crop}") + + return slice(start, stop, step) + + +def main(): + parser = argparse.ArgumentParser(description="Crop images to a region of interest", + epilog="If a start or stop value is omitted, the start or end of the axis is assumed.") + parser.add_argument("-i", metavar="I1:I2[:-1]", + help="Start/stop [flip] along first axis (0-indexed)") + parser.add_argument("-j", metavar="J1:J2[:-1]", + help="Start/stop [flip] along second axis (0-indexed)") + parser.add_argument("-k", metavar="K1:K2[:-1]", + help="Start/stop [flip] along third axis (0-indexed)") + parser.add_argument("-t", metavar="T1:T2", help="Start/stop along fourth axis (0-indexed)") + parser.add_argument("in_file", help="Image file to crop") + parser.add_argument("out_file", help="Output file name") + + opts = parser.parse_args() + + try: + islice = parse_slice(opts.i) + jslice = parse_slice(opts.j) + kslice = parse_slice(opts.k) + tslice = parse_slice(opts.t, allow_step=False) + except ValueError as err: + print(f"Could not parse input arguments. Reason follows.\n{err}") + return 1 + + img = nb.load(opts.in_file) + try: + sliced_img = lossless_slice(img, (islice, jslice, kslice, tslice)[:img.ndim]) + except: + print("Could not slice image. Full traceback follows.") + raise + nb.save(sliced_img, opts.out_file) + return 0 diff --git a/setup.cfg b/setup.cfg index 1231a7502e..287ae64a04 100644 --- a/setup.cfg +++ b/setup.cfg @@ -77,6 +77,7 @@ console_scripts = nib-nifti-dx=nibabel.cmdline.nifti_dx:main nib-tck2trk=nibabel.cmdline.tck2trk:main nib-trk2tck=nibabel.cmdline.trk2tck:main + nib-roi=nibabel.cmdline.roi:main parrec2nii=nibabel.cmdline.parrec2nii:main [options.package_data] From ab960440f9b41219b10850651f5e64d2a4fefab6 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 15 Sep 2020 09:08:37 -0400 Subject: [PATCH 2/7] ENH: Cover edge cases like non-scaling images, bad slices --- nibabel/cmdline/roi.py | 50 +++++++++++++++++++++++++++++++----------- 1 file changed, 37 insertions(+), 13 deletions(-) diff --git a/nibabel/cmdline/roi.py b/nibabel/cmdline/roi.py index 485548c08e..ad8305af93 100644 --- a/nibabel/cmdline/roi.py +++ b/nibabel/cmdline/roi.py @@ -1,3 +1,5 @@ +import sys +import os import argparse import nibabel as nb @@ -6,11 +8,13 @@ def lossless_slice(img, slicers): if not nb.imageclasses.spatial_axes_first(img): raise ValueError("Cannot slice an image that is not known to have spatial axes first") - roi_img = img.__class__( - img.dataobj._get_unscaled(slicers), - affine=img.slicer.slice_affine(slicers), - header=img.header) - roi_img.header.set_slope_inter(img.dataobj.slope, img.dataobj.inter) + scaling = hasattr(img.dataobj, 'slope') + + data = img.dataobj._get_unscaled(slicers) if scaling else img.dataobj[slicers] + roi_img = img.__class__(data, affine=img.slicer.slice_affine(slicers), header=img.header) + + if scaling: + roi_img.header.set_slope_inter(img.dataobj.slope, img.dataobj.inter) return roi_img @@ -20,7 +24,7 @@ def parse_slice(crop, allow_step=True): start, stop, *extra = [int(val) if val else None for val in crop.split(":")] if len(extra) > 1: raise ValueError(f"Cannot parse specification: {crop}") - if extra and not allow_step: + if not allow_step and extra and extra[0] not in (1, None): raise ValueError(f"Step entry not permitted: {crop}") step = extra[0] if extra else None @@ -30,9 +34,19 @@ def parse_slice(crop, allow_step=True): return slice(start, stop, step) -def main(): - parser = argparse.ArgumentParser(description="Crop images to a region of interest", - epilog="If a start or stop value is omitted, the start or end of the axis is assumed.") +def sanitize(args): + # Argparse likes to treat "-1:..." as a flag + return [f' {arg}' if arg[0] == '-' and ":" in arg else arg + for arg in args] + + +def main(args=None): + if args is None: + args = sys.argv[1:] + parser = argparse.ArgumentParser( + description="Crop images to a region of interest", + epilog="If a start or stop value is omitted, the start or end of the axis is assumed.") + parser.add_argument('--version', action='version', version=nb.__version__) parser.add_argument("-i", metavar="I1:I2[:-1]", help="Start/stop [flip] along first axis (0-indexed)") parser.add_argument("-j", metavar="J1:J2[:-1]", @@ -43,7 +57,7 @@ def main(): parser.add_argument("in_file", help="Image file to crop") parser.add_argument("out_file", help="Output file name") - opts = parser.parse_args() + opts = parser.parse_args(args=sanitize(args)) try: islice = parse_slice(opts.i) @@ -54,10 +68,20 @@ def main(): print(f"Could not parse input arguments. Reason follows.\n{err}") return 1 - img = nb.load(opts.in_file) + kwargs = {} + if os.path.realpath(opts.in_file) == os.path.realpath(opts.out_file): + kwargs['mmap'] = False + img = nb.load(opts.in_file, **kwargs) + + slicers = (islice, jslice, kslice, tslice)[:img.ndim] + expected_shape = nb.fileslice.predict_shape(slicers, img.shape) + if any(dim == 0 for dim in expected_shape): + print(f"Cannot take zero-length slices. Predicted shape {expected_shape}.") + return 1 + try: - sliced_img = lossless_slice(img, (islice, jslice, kslice, tslice)[:img.ndim]) - except: + sliced_img = lossless_slice(img, slicers) + except Exception: print("Could not slice image. Full traceback follows.") raise nb.save(sliced_img, opts.out_file) From 310a242539c5ae45d1364ada5cdca564a0c38be6 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 15 Sep 2020 09:09:00 -0400 Subject: [PATCH 3/7] TEST: Test nib-roi and helpers --- nibabel/cmdline/tests/test_roi.py | 115 ++++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) create mode 100644 nibabel/cmdline/tests/test_roi.py diff --git a/nibabel/cmdline/tests/test_roi.py b/nibabel/cmdline/tests/test_roi.py new file mode 100644 index 0000000000..ddd105a30c --- /dev/null +++ b/nibabel/cmdline/tests/test_roi.py @@ -0,0 +1,115 @@ +import os +import numpy as np +import nibabel as nb +from nibabel.cmdline.roi import lossless_slice, parse_slice, main +from nibabel.testing import data_path + +import unittest +import pytest + + +def test_parse_slice(): + assert parse_slice(None) == slice(None) + assert parse_slice("1:5") == slice(1, 5) + assert parse_slice("1:") == slice(1, None) + assert parse_slice(":5") == slice(None, 5) + assert parse_slice(":-1") == slice(None, -1) + assert parse_slice("-5:-1") == slice(-5, -1) + assert parse_slice("1:5:") == slice(1, 5, None) + assert parse_slice("1::") == slice(1, None, None) + assert parse_slice(":5:") == slice(None, 5, None) + assert parse_slice(":-1:") == slice(None, -1, None) + assert parse_slice("-5:-1:") == slice(-5, -1, None) + assert parse_slice("1:5:1") == slice(1, 5, 1) + assert parse_slice("1::1") == slice(1, None, 1) + assert parse_slice(":5:1") == slice(None, 5, 1) + assert parse_slice(":-1:1") == slice(None, -1, 1) + assert parse_slice("-5:-1:1") == slice(-5, -1, 1) + assert parse_slice("5:1:-1") == slice(5, 1, -1) + assert parse_slice(":1:-1") == slice(None, 1, -1) + assert parse_slice("5::-1") == slice(5, None, -1) + assert parse_slice("-1::-1") == slice(-1, None, -1) + assert parse_slice("-1:-5:-1") == slice(-1, -5, -1) + + # Max of start:stop:step + with pytest.raises(ValueError): + parse_slice("1:2:3:4") + # Integers only + with pytest.raises(ValueError): + parse_slice("abc:2:3") + with pytest.raises(ValueError): + parse_slice("1.2:2:3") + # Unit steps only + with pytest.raises(ValueError): + parse_slice("1:5:2") + + +def test_parse_slice_disallow_step(): + # Permit steps of 1 + assert parse_slice("1:5", False) == slice(1, 5) + assert parse_slice("1:5:", False) == slice(1, 5) + assert parse_slice("1:5:1", False) == slice(1, 5, 1) + # Disable other steps + with pytest.raises(ValueError): + parse_slice("1:5:-1", False) + with pytest.raises(ValueError): + parse_slice("1:5:-2", False) + + +def test_lossless_slice_unknown_axes(): + img = nb.load(os.path.join(data_path, 'minc1_4d.mnc')) + with pytest.raises(ValueError): + lossless_slice(img, (slice(None), slice(None), slice(None))) + + +def test_lossless_slice_scaling(tmp_path): + fname = tmp_path / 'image.nii' + img = nb.Nifti1Image(np.random.uniform(-20000, 20000, (5, 5, 5, 5)), affine=np.eye(4)) + img.header.set_data_dtype("int16") + img.to_filename(fname) + img1 = nb.load(fname) + sliced_fname = tmp_path / 'sliced.nii' + lossless_slice(img1, (slice(None), slice(None), slice(2, 4))).to_filename(sliced_fname) + img2 = nb.load(sliced_fname) + + assert np.array_equal(img1.get_fdata()[:, :, 2:4], img2.get_fdata()) + assert np.array_equal(img1.dataobj.get_unscaled()[:, :, 2:4], img2.dataobj.get_unscaled()) + assert img1.dataobj.slope == img2.dataobj.slope + assert img1.dataobj.inter == img2.dataobj.inter + + +@pytest.mark.parametrize("inplace", (True, False)) +def test_nib_roi(tmp_path, inplace): + in_file = os.path.join(data_path, 'functional.nii') + out_file = str(tmp_path / 'sliced.nii') + in_img = nb.load(in_file) + + if inplace: + in_img.to_filename(out_file) + in_file = out_file + + retval = main([in_file, out_file, '-i', '1:-1', '-j', '-1:1:-1', '-k', '::', '-t', ':5']) + assert retval == 0 + + out_img = nb.load(out_file) + in_data = in_img.dataobj[:] + in_sliced = in_img.slicer[1:-1, -1:1:-1, :, :5] + assert out_img.shape == in_sliced.shape + assert np.array_equal(in_data[1:-1, -1:1:-1, :, :5], out_img.dataobj) + assert np.allclose(in_sliced.dataobj, out_img.dataobj) + assert np.allclose(in_sliced.affine, out_img.affine) + + +@pytest.mark.parametrize("args, errmsg", ( + (("-i", "1:1"), "Cannot take zero-length slice"), + (("-j", "1::2"), "Downsampling is not supported"), + (("-t", "5::-1"), "Step entry not permitted"), +)) +def test_nib_roi_bad_slices(capsys, args, errmsg): + in_file = os.path.join(data_path, 'functional.nii') + + retval = main([in_file, '/dev/null', *args]) + assert retval != 0 + captured = capsys.readouterr() + assert errmsg in captured.out + From e969374d30dbfcc64fb88848dc1ec6467fa1a517 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 23 Sep 2020 09:11:00 -0400 Subject: [PATCH 4/7] Update nibabel/cmdline/tests/test_roi.py --- nibabel/cmdline/tests/test_roi.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/nibabel/cmdline/tests/test_roi.py b/nibabel/cmdline/tests/test_roi.py index ddd105a30c..46c08113c6 100644 --- a/nibabel/cmdline/tests/test_roi.py +++ b/nibabel/cmdline/tests/test_roi.py @@ -108,8 +108,7 @@ def test_nib_roi(tmp_path, inplace): def test_nib_roi_bad_slices(capsys, args, errmsg): in_file = os.path.join(data_path, 'functional.nii') - retval = main([in_file, '/dev/null', *args]) + retval = main([in_file, os.devnull, *args]) assert retval != 0 captured = capsys.readouterr() assert errmsg in captured.out - From d0a64254924a4a56119369bdb7a3ff678fd6603c Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 23 Sep 2020 10:42:29 -0400 Subject: [PATCH 5/7] TEST: Check MINC-2, Analyze, and main() --- nibabel/cmdline/tests/test_roi.py | 39 +++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/nibabel/cmdline/tests/test_roi.py b/nibabel/cmdline/tests/test_roi.py index 46c08113c6..6e349d65db 100644 --- a/nibabel/cmdline/tests/test_roi.py +++ b/nibabel/cmdline/tests/test_roi.py @@ -5,6 +5,7 @@ from nibabel.testing import data_path import unittest +from unittest import mock import pytest @@ -78,6 +79,23 @@ def test_lossless_slice_scaling(tmp_path): assert img1.dataobj.inter == img2.dataobj.inter +def test_lossless_slice_noscaling(tmp_path): + fname = tmp_path / 'image.img' + img = nb.AnalyzeImage(np.random.uniform(-20000, 20000, (5, 5, 5, 5)).astype("float32"), + affine=np.eye(4)) + img.header.set_data_dtype("float32") + img.to_filename(fname) + img1 = nb.load(fname) + sliced_fname = tmp_path / 'sliced.img' + lossless_slice(img1, (slice(None), slice(None), slice(2, 4))).to_filename(sliced_fname) + img2 = nb.load(sliced_fname) + + assert np.array_equal(img1.get_fdata()[:, :, 2:4], img2.get_fdata()) + assert np.array_equal(img1.dataobj.get_unscaled()[:, :, 2:4], img2.dataobj.get_unscaled()) + assert img1.dataobj.slope == img2.dataobj.slope + assert img1.dataobj.inter == img2.dataobj.inter + + @pytest.mark.parametrize("inplace", (True, False)) def test_nib_roi(tmp_path, inplace): in_file = os.path.join(data_path, 'functional.nii') @@ -112,3 +130,24 @@ def test_nib_roi_bad_slices(capsys, args, errmsg): assert retval != 0 captured = capsys.readouterr() assert errmsg in captured.out + + +def test_entrypoint(capsys): + # Check that we handle missing args as expected + with mock.patch("sys.argv", ["nib-roi", "--help"]): + try: + retval = main() + except SystemExit: + pass + else: + assert False, "argparse exits on --help. If changing to another parser, update test." + captured = capsys.readouterr() + assert captured.out.startswith("usage: nib-roi") + + +def test_nib_roi_unknown_axes(capsys): + in_file = os.path.join(data_path, 'minc1_4d.mnc') + with pytest.raises(ValueError): + main([in_file, os.devnull, "-i", ":"]) + captured = capsys.readouterr() + assert "Could not slice image." in captured.out From 4e985a7969a7e79d1c3d6d5b01ceb07cf25d58fd Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Fri, 9 Oct 2020 08:39:50 -0400 Subject: [PATCH 6/7] TEST: Use MGHImage to avoid scipy dependency --- nibabel/cmdline/tests/test_roi.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/nibabel/cmdline/tests/test_roi.py b/nibabel/cmdline/tests/test_roi.py index 6e349d65db..4c640e9136 100644 --- a/nibabel/cmdline/tests/test_roi.py +++ b/nibabel/cmdline/tests/test_roi.py @@ -80,13 +80,12 @@ def test_lossless_slice_scaling(tmp_path): def test_lossless_slice_noscaling(tmp_path): - fname = tmp_path / 'image.img' - img = nb.AnalyzeImage(np.random.uniform(-20000, 20000, (5, 5, 5, 5)).astype("float32"), - affine=np.eye(4)) - img.header.set_data_dtype("float32") + fname = tmp_path / 'image.mgh' + img = nb.MGHImage(np.random.uniform(-20000, 20000, (5, 5, 5, 5)).astype("float32"), + affine=np.eye(4)) img.to_filename(fname) img1 = nb.load(fname) - sliced_fname = tmp_path / 'sliced.img' + sliced_fname = tmp_path / 'sliced.mgh' lossless_slice(img1, (slice(None), slice(None), slice(2, 4))).to_filename(sliced_fname) img2 = nb.load(sliced_fname) From 4623e12d681d2127ec5362268ab3301eed88020d Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Fri, 9 Oct 2020 08:40:03 -0400 Subject: [PATCH 7/7] FIX: Improved test for scaling --- nibabel/cmdline/roi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/cmdline/roi.py b/nibabel/cmdline/roi.py index ad8305af93..0631ecc0d1 100644 --- a/nibabel/cmdline/roi.py +++ b/nibabel/cmdline/roi.py @@ -8,7 +8,7 @@ def lossless_slice(img, slicers): if not nb.imageclasses.spatial_axes_first(img): raise ValueError("Cannot slice an image that is not known to have spatial axes first") - scaling = hasattr(img.dataobj, 'slope') + scaling = hasattr(img.header, 'set_slope_inter') data = img.dataobj._get_unscaled(slicers) if scaling else img.dataobj[slicers] roi_img = img.__class__(data, affine=img.slicer.slice_affine(slicers), header=img.header)