Skip to content
Open
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .python-version
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
3.12
>=3.10, < 3.13
>=3.11, < 3.13
6 changes: 3 additions & 3 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,11 @@ build-backend = "setuptools.build_meta"

[project]
name = "scilpy"
version = "2.2.2"
version = "2.3.0"
description = "Scilpy: diffusion MRI tools and utilities"
authors = [{ name = "SCIL Team" }]
readme = "README.md"
requires-python = ">=3.10, <3.13"
requires-python = ">=3.11, <3.13"
license-files = ["LICENSE"]
classifiers = [
"Development Status :: 3 - Alpha",
Expand All @@ -31,7 +31,7 @@ dependencies = [
"cvxpy==1.6.*",
"cycler==0.12.*",
"deepdiff>=8.1,<8.7",
"dipy==1.11.*",
"dipy==1.12.*",
"dmri-amico==2.1.*",
"dmri-commit>=2.4.2,<2.5.0",
"docopt==0.6.*",
Expand Down
25 changes: 18 additions & 7 deletions src/scilpy/cli/scil_denoising_nlmeans.py
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,10 @@ def main():
if args.sigma is not None:
logging.info('User supplied noise standard deviation is {}'
.format(args.sigma))
sigma = np.ones(vol_data.shape[:3]) * args.sigma
if nb_volumes > 1:
sigma = np.full(nb_volumes, args.sigma, dtype=np.float32)
else:
sigma = float(args.sigma)
elif args.basic_sigma:
if args.mask_sigma:
mask_sigma = get_data_as_mask(nib.load(args.mask_sigma))
Expand All @@ -206,8 +209,8 @@ def main():
sigma = np.median(sigma) # Managing 4D data.
logging.info('The median noise is: {}'.format(sigma))

# Broadcast the single value to a whole 3D volume for nlmeans
sigma = np.ones(vol_data.shape) * sigma
if nb_volumes > 1:
sigma = np.full(nb_volumes, sigma, dtype=np.float32)
else: # --piesno
logging.info("Computing sigma: one value per slice.")
sigma, mask_noise = estimate_piesno_sigma(vol_data, args.number_coils)
Expand All @@ -219,14 +222,22 @@ def main():
header=vol.header),
args.save_piesno_mask)

# Broadcast the values per slice to a whole 3D volume for nlmeans
# Keep a 3D sigma map (one value per slice) for PIESNO.
sigma = np.ones(vol_data.shape[:3]) * sigma[None, None, :]

with warnings.catch_warnings():
warnings.simplefilter("ignore", category=DeprecationWarning)
data_denoised = nlmeans(
vol_data, sigma, mask=mask_denoise, rician=not args.gaussian,
num_threads=args.nbr_processes)
if vol_data.ndim == 4 and np.ndim(sigma) == 3:
data_denoised = np.empty_like(vol_data)
for vol_idx in range(vol_data.shape[-1]):
data_denoised[..., vol_idx] = nlmeans(
vol_data[..., vol_idx], sigma,
mask=mask_denoise, rician=not args.gaussian,
num_threads=args.nbr_processes)
else:
data_denoised = nlmeans(
vol_data, sigma, mask=mask_denoise, rician=not args.gaussian,
num_threads=args.nbr_processes)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this change would be because nlmeans does not support 4D volumes anymore?


# Saving
nib.save(nib.Nifti1Image(data_denoised, vol.affine, header=vol.header),
Expand Down
8 changes: 4 additions & 4 deletions src/scilpy/cli/tests/test_tractogram_convert.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,10 @@ def test_help_option(script_runner):

def test_execution_surface_vtk_fib(script_runner, monkeypatch):
monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
in_fib = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
'gyri_fanning.fib')
in_trk = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
'gyri_fanning.trk')
in_fa = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
'fa.nii.gz')
ret = script_runner.run(['scil_tractogram_convert', in_fib,
'gyri_fanning.trk', '--reference', in_fa])
ret = script_runner.run(['scil_tractogram_convert', in_trk,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I take it we can't convert fib to trk anymore?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can but the fib file we have for our test are not correct.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we update the file then, before merging the PR?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@frheault do you think it's worth it ? we can remove the file or I create a new/fixed .fib file by converting the trk to fib.

'gyri_fanning converted.fib', '--reference', in_fa])
assert ret.success
8 changes: 4 additions & 4 deletions src/scilpy/cli/tests/test_tractogram_flip.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,10 @@ def test_help_option(script_runner):

def test_execution_surface_vtk_fib(script_runner, monkeypatch):
monkeypatch.chdir(os.path.expanduser(tmp_dir.name))
in_fib = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
'gyri_fanning.fib')
in_trk = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not sure the name of this test really fits the content? What is surface_vtk_fib?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"surface_vtk_fib" is the test folder used here.

'gyri_fanning.trk')
in_fa = os.path.join(SCILPY_HOME, 'surface_vtk_fib',
'fa.nii.gz')
ret = script_runner.run(['scil_tractogram_flip', in_fib,
'gyri_fanning.tck', 'x', '--reference', in_fa])
ret = script_runner.run(['scil_tractogram_flip', in_trk,
'gyri_fanning_fliped.tck', 'x', '--reference', in_fa])
assert ret.success
15 changes: 15 additions & 0 deletions src/scilpy/tractanalysis/tests/test_voxel_boundary_intersection.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,3 +38,18 @@ def test_subdivide():

# streamline 3 should be split into 11 (12 points)
assert split[2].shape[0] == 12


def test_subdivide_float64_streamlines():
streamline = np.asarray([[0.0, 0.0, 0.0],
[1.2, 0.5, 0.0]], dtype=np.float64)

fake_ref = nib.Nifti1Image(np.zeros((3, 3, 3)), affine=np.eye(4))
sft = StatefulTractogram([streamline], reference=fake_ref,
origin=Origin('corner'), space=Space('vox'))

split = subdivide_streamlines_at_voxel_faces(sft.streamlines)

assert split[0].dtype == np.float32
assert np.allclose(split[0][0], streamline[0])
assert np.allclose(split[0][-1], streamline[-1])
7 changes: 5 additions & 2 deletions src/scilpy/tractanalysis/voxel_boundary_intersection.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -45,13 +45,16 @@ def subdivide_streamlines_at_voxel_faces(streamlines):
Updated streamline coordinates with added coordinate points
at voxel boundaries.
"""
streamlines_data = np.ascontiguousarray(streamlines._data,
dtype=np.float32)

cdef:
cnp.npy_intp nb_streamlines = len(streamlines._lengths)
cnp.npy_intp at_point = 0

# Multiplying by 6 is simply a heuristic to avoiding resizing too many
# times. In my bundles tests, I had either 0 or 1 resize.
cnp.npy_intp max_points = (streamlines.get_data().size // 6) * 12
cnp.npy_intp max_points = (streamlines_data.size // 6) * 12

new_array_sequence = nib.streamlines.array_sequence.ArraySequence()
new_array_sequence._lengths.resize(nb_streamlines)
Expand All @@ -61,7 +64,7 @@ def subdivide_streamlines_at_voxel_faces(streamlines):
cdef:
cnp.npy_intp[:] lengths_view_in = streamlines._lengths
cnp.npy_intp[:] offsets_view_in = streamlines._offsets
float[:, :] data_view_in = streamlines._data
float[:, :] data_view_in = streamlines_data
cnp.npy_intp[:] lengths_view_out = new_array_sequence._lengths
cnp.npy_intp[:] offsets_view_out = new_array_sequence._offsets
cnp.float32_t[:] data_view_out = new_array_sequence._data
Expand Down
Loading