Skip to content

Commit

Permalink
Merge pull request #1020 from frheault/lesions_divide
Browse files Browse the repository at this point in the history
Modify lesions_info to support labels instead of mask
  • Loading branch information
arnaudbore authored Aug 7, 2024
2 parents 2d8520c + ca306f8 commit 7b99033
Show file tree
Hide file tree
Showing 4 changed files with 37 additions and 16 deletions.
19 changes: 18 additions & 1 deletion scilpy/image/labels.py
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,8 @@ def get_binary_mask_from_labels(atlas, label_list):
return mask


def get_labels_from_mask(mask_data, labels=None, background_label=0):
def get_labels_from_mask(mask_data, labels=None, background_label=0,
min_voxel_count=0):
"""
Get labels from a binary mask which contains multiple blobs. Each blob
will be assigned a label, by default starting from 1. Background will
Expand All @@ -83,6 +84,9 @@ def get_labels_from_mask(mask_data, labels=None, background_label=0):
label.
background_label: int
Label for the background.
min_voxel_count: int, optional
Minimum number of voxels for a blob to be considered. Blobs with fewer
voxels will be ignored.
Returns
-------
Expand All @@ -91,6 +95,19 @@ def get_labels_from_mask(mask_data, labels=None, background_label=0):
"""
# Get the number of structures and assign labels to each blob
label_map, nb_structures = ndi.label(mask_data)
if min_voxel_count:
new_count = 0
for label in range(1, nb_structures + 1):
if np.count_nonzero(label_map == label) < min_voxel_count:
label_map[label_map == label] = 0
else:
new_count += 1
label_map[label_map == label] = new_count
logging.debug(
f"Ignored blob {nb_structures-new_count} with fewer "
"than {min_voxel_count} voxels")
nb_structures = new_count

# Assign labels to each blob if provided
if labels:
# Only keep the first nb_structures labels if the number of labels
Expand Down
10 changes: 9 additions & 1 deletion scripts/scil_labels_from_mask.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,10 @@ def _build_arg_parser():
p.add_argument('--background_label', default=0, type=int,
help='Label to assign to the background. [%(default)s]')

p.add_argument('--min_volume', type=float, default=7,
help='Minimum volume in mm3 [%(default)s],'
'Useful for lesions.')

add_verbose_arg(p)
add_overwrite_arg(p)

Expand All @@ -50,9 +54,13 @@ def main():
# Load mask and get data
mask_img = nib.load(args.in_mask)
mask_data = get_data_as_mask(mask_img)
voxel_volume = np.prod(np.diag(mask_img.affine)[:3])
min_voxel_count = args.min_volume // voxel_volume

# Get labels from mask
label_map = get_labels_from_mask(
mask_data, args.labels, args.background_label)
mask_data, args.labels, args.background_label,
min_voxel_count=min_voxel_count)
# Save result
out_img = nib.Nifti1Image(label_map.astype(np.uint16), mask_img.affine)
nib.save(out_img, args.out_labels)
Expand Down
22 changes: 9 additions & 13 deletions scripts/scil_lesions_info.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,12 @@

"""
This script will output informations about lesion load in bundle(s).
The input can either be streamlines, binary bundle map, or a bundle voxel
label map.
The input lesion file is a labeled volume (.nii.gz) where each lesion is
represented by a unique label. To label a lesion file use
scil_labels_from_mask.py
To be considered a valid lesion, the lesion volume must be at least
min_lesion_vol mm3. This avoid the detection of thousand of single voxel
lesions if an automatic lesion segmentation tool is used.
Then, the second input can either be streamlines, binary bundle mask, or a
bundle voxel label map.
Formerly: scil_analyse_lesions_load.py
"""
Expand All @@ -20,7 +20,6 @@

import nibabel as nib
import numpy as np
import scipy.ndimage as ndi

from scilpy.image.labels import get_data_as_labels
from scilpy.io.image import get_data_as_mask
Expand All @@ -40,7 +39,7 @@ def _build_arg_parser():
formatter_class=argparse.RawTextHelpFormatter)

p.add_argument('in_lesion',
help='Binary mask of the lesion(s) (.nii.gz).')
help='Lesions file as labels (.nii.gz).')
p.add_argument('out_json',
help='Output file for lesion information (.json).')
p1 = p.add_mutually_exclusive_group()
Expand All @@ -51,8 +50,6 @@ def _build_arg_parser():
p1.add_argument('--bundle_labels_map',
help='Path of the bundle labels map (.nii.gz).')

p.add_argument('--min_lesion_vol', type=float, default=7,
help='Minimum lesion volume in mm3 [%(default)s].')
p.add_argument('--out_lesion_atlas', metavar='FILE',
help='Save the labelized lesion(s) map (.nii.gz).')
p.add_argument('--out_lesion_stats', metavar='FILE',
Expand Down Expand Up @@ -89,7 +86,8 @@ def main():
reference=args.reference)

lesion_img = nib.load(args.in_lesion)
lesion_data = get_data_as_mask(lesion_img, dtype=bool)
voxel_sizes = lesion_img.header.get_zooms()[0:3]
lesion_atlas = get_data_as_labels(lesion_img)

if args.bundle:
bundle_name, _ = split_name_with_nii(os.path.basename(args.bundle))
Expand All @@ -98,7 +96,7 @@ def main():
sft.to_corner()
streamlines = sft.get_streamlines_copy()
map_data = compute_tract_counts_map(streamlines,
lesion_data.shape)
lesion_atlas.shape)
map_data[map_data > 0] = 1
elif args.bundle_mask:
bundle_name, _ = split_name_with_nii(
Expand All @@ -112,8 +110,6 @@ def main():
map_data = get_data_as_labels(map_img)

is_single_label = args.bundle_labels_map is None
voxel_sizes = lesion_img.header.get_zooms()[0:3]
lesion_atlas, _ = ndi.label(lesion_data)

lesion_load_dict = compute_lesion_stats(
map_data, lesion_atlas, single_label=is_single_label,
Expand Down
2 changes: 1 addition & 1 deletion scripts/tests/test_labels_from_mask.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ def test_execution(script_runner, monkeypatch):
'bundle_4_head_tail_offset.nii.gz')
ret = script_runner.run('scil_labels_from_mask.py',
in_mask, 'labels_from_mask.nii.gz',
'-f')
'--min_volume', '0', '-f')
assert ret.success


Expand Down

0 comments on commit 7b99033

Please sign in to comment.