From 29caff4024110884c5846dababa0c0ad4133cc5f Mon Sep 17 00:00:00 2001 From: George Aidins Date: Thu, 31 Aug 2023 15:08:58 -0400 Subject: [PATCH 1/3] Large FOV bugfix, add pre/post steps to crop image to ICV mask before applying DLMUSE --- niCHARTPipelines/CombineMasks.py | 50 ++++++++++++++++++++++++++++---- niCHARTPipelines/MaskImage.py | 43 +++++++++++++++++++++++++-- 2 files changed, 86 insertions(+), 7 deletions(-) diff --git a/niCHARTPipelines/CombineMasks.py b/niCHARTPipelines/CombineMasks.py index 22617da..9d19c87 100644 --- a/niCHARTPipelines/CombineMasks.py +++ b/niCHARTPipelines/CombineMasks.py @@ -1,7 +1,36 @@ +import numpy as np import nibabel as nib -from nibabel.orientations import axcodes2ornt, inv_ornt_aff, ornt_transform from nipype.interfaces.image import Reorient +from nibabel.orientations import axcodes2ornt, ornt_transform, inv_ornt_aff +## Find bounding box for the foreground values in img, with a given padding percentage +def calc_bbox_with_padding(img, perc_pad = 10): + + ## Output is the coordinates of the bounding box + bcoors = np.zeros([3,2], dtype=int) + + ## Find coors in each axis + for sel_axis in [0, 1, 2]: + + ## Get axes other than the selected + other_axes = [0, 1, 2] + other_axes.remove(sel_axis) + + ## Get img dim in selected axis + dim = img.shape[sel_axis] + + ## Find bounding box (index of first and last non-zero slices) + nonzero = np.any(img, axis = tuple(other_axes)) + bbox= np.where(nonzero)[0][[0,-1]] + + ## Add padding + size_pad = int(np.round((bbox[1] - bbox[0]) * perc_pad / 100)) + b_min = int(np.max([0, bbox[0] - size_pad])) + b_max = int(np.min([dim, bbox[1] + size_pad])) + + bcoors[sel_axis, :] = [b_min, b_max] + + return bcoors def apply_combine(in_img_name, icv_img_name, out_img_name): '''Combine icv and muse masks. @@ -12,13 +41,24 @@ def apply_combine(in_img_name, icv_img_name, out_img_name): img_in = nii_in.get_fdata() img_icv = nii_icv.get_fdata() + + ################################ + ## INFO: nnunet hallucinated on images with large FOV. To solve this problem + ## we added pre/post processing steps to crop initial image around ICV + ## mask before sending to DLMUSE + ## + ## MUSE image (img_in) may be cropped. Pad it to initial image size + bcoors = calc_bbox_with_padding(img_icv) + img_out = img_icv * 0 + img_out[bcoors[0,0]:bcoors[0,1], bcoors[1,0]:bcoors[1,1], bcoors[2,0]:bcoors[2,1]] = img_in + ################################ # Merge masks : Add a new label (1) to MUSE for foreground voxels in ICV that is not in MUSE # this label will mainly represent cortical CSF - img_in[(img_in==0) & (img_icv>0)] = 1 + img_out[(img_out==0) & (img_icv>0)] = 1 - img_in = img_in.astype(int) + img_out = img_out.astype(int) ## Save out image - nii_out = nib.Nifti1Image(img_in, nii_in.affine, nii_in.header) - nii_out.to_filename(out_img_name) + nii_out = nib.Nifti1Image(img_out, nii_in.affine, nii_in.header) + nii_out.to_filename(out_img_name) \ No newline at end of file diff --git a/niCHARTPipelines/MaskImage.py b/niCHARTPipelines/MaskImage.py index 819f513..64c9eea 100644 --- a/niCHARTPipelines/MaskImage.py +++ b/niCHARTPipelines/MaskImage.py @@ -1,5 +1,35 @@ +import numpy as np import nibabel as nib +## Find bounding box for the foreground values in img, with a given padding percentage +def calc_bbox_with_padding(img, perc_pad = 10): + + ## Output is the coordinates of the bounding box + bcoors = np.zeros([3,2], dtype=int) + + ## Find coors in each axis + for sel_axis in [0, 1, 2]: + + ## Get axes other than the selected + other_axes = [0, 1, 2] + other_axes.remove(sel_axis) + + ## Get img dim in selected axis + dim = img.shape[sel_axis] + + ## Find bounding box (index of first and last non-zero slices) + nonzero = np.any(img, axis = tuple(other_axes)) + bbox= np.where(nonzero)[0][[0,-1]] + + ## Add padding + size_pad = int(np.round((bbox[1] - bbox[0]) * perc_pad / 100)) + b_min = int(np.max([0, bbox[0] - size_pad])) + b_max = int(np.min([dim, bbox[1] + size_pad])) + + bcoors[sel_axis, :] = [b_min, b_max] + + return bcoors + ###---------mask image----------- def apply_mask(in_img_name, mask_img_name, out_img_name): @@ -13,6 +43,15 @@ def apply_mask(in_img_name, mask_img_name, out_img_name): ## Mask image img_in[img_mask == 0] = 0 + ################################ + ## INFO: nnunet hallucinated on images with large FOV. To solve this problem + ## we added pre/post processing steps to crop initial image around ICV + ## mask before sending to DLMUSE + ## + ## Crop image + bcoors = calc_bbox_with_padding(img_mask) + img_in_crop = img_in[bcoors[0,0]:bcoors[0,1], bcoors[1,0]:bcoors[1,1], bcoors[2,0]:bcoors[2,1]] + ## Save out image - nii_out = nib.Nifti1Image(img_in, nii_in.affine, nii_in.header) - nii_out.to_filename(out_img_name) + nii_out = nib.Nifti1Image(img_in_crop, nii_in.affine, nii_in.header) + nii_out.to_filename(out_img_name) \ No newline at end of file From db1c59e045f90bf5b93954a58b885dcf52588030 Mon Sep 17 00:00:00 2001 From: George Aidinis Date: Fri, 1 Sep 2023 14:28:14 -0400 Subject: [PATCH 2/3] Bugfix for bounding box estimation --- niCHARTPipelines/CombineMasks.py | 24 ++++++++++++++++++++---- niCHARTPipelines/MaskImage.py | 19 +++++++++++++++++-- 2 files changed, 37 insertions(+), 6 deletions(-) diff --git a/niCHARTPipelines/CombineMasks.py b/niCHARTPipelines/CombineMasks.py index 9d19c87..ccf0f30 100644 --- a/niCHARTPipelines/CombineMasks.py +++ b/niCHARTPipelines/CombineMasks.py @@ -1,14 +1,30 @@ -import numpy as np import nibabel as nib +import numpy as np +from nibabel.orientations import axcodes2ornt, inv_ornt_aff, ornt_transform from nipype.interfaces.image import Reorient -from nibabel.orientations import axcodes2ornt, ornt_transform, inv_ornt_aff +from scipy import ndimage +from scipy.ndimage.measurements import label + ## Find bounding box for the foreground values in img, with a given padding percentage def calc_bbox_with_padding(img, perc_pad = 10): + img = img.astype('uint8') + ## Output is the coordinates of the bounding box bcoors = np.zeros([3,2], dtype=int) + ## Find the largest connected component + ## INFO: In images with very large FOV DLICV may have small isolated regions in + ## boundaries; so we calculate the bounding box based on the brain, not all + ## foreground voxels + str_3D = np.array([[[0, 0, 0], [0, 1, 0], [0, 0, 0]], + [[0, 1, 0], [1, 1, 1], [0, 1, 0]], + [[0, 0, 0], [0, 1, 0], [0, 0, 0]]], dtype='uint8') + labeled, ncomp = label(img, str_3D) + sizes = ndimage.sum(img, labeled, range(ncomp + 1)) + img_largest_cc = (labeled == np.argmax(sizes)).astype(int) + ## Find coors in each axis for sel_axis in [0, 1, 2]: @@ -17,10 +33,10 @@ def calc_bbox_with_padding(img, perc_pad = 10): other_axes.remove(sel_axis) ## Get img dim in selected axis - dim = img.shape[sel_axis] + dim = img_largest_cc.shape[sel_axis] ## Find bounding box (index of first and last non-zero slices) - nonzero = np.any(img, axis = tuple(other_axes)) + nonzero = np.any(img_largest_cc, axis = tuple(other_axes)) bbox= np.where(nonzero)[0][[0,-1]] ## Add padding diff --git a/niCHARTPipelines/MaskImage.py b/niCHARTPipelines/MaskImage.py index 64c9eea..6f054a5 100644 --- a/niCHARTPipelines/MaskImage.py +++ b/niCHARTPipelines/MaskImage.py @@ -1,12 +1,27 @@ import numpy as np import nibabel as nib +from scipy import ndimage +from scipy.ndimage.measurements import label ## Find bounding box for the foreground values in img, with a given padding percentage def calc_bbox_with_padding(img, perc_pad = 10): + img = img.astype('uint8') + ## Output is the coordinates of the bounding box bcoors = np.zeros([3,2], dtype=int) + ## Find the largest connected component + ## INFO: In images with very large FOV DLICV may have small isolated regions in + ## boundaries; so we calculate the bounding box based on the brain, not all + ## foreground voxels + str_3D = np.array([[[0, 0, 0], [0, 1, 0], [0, 0, 0]], + [[0, 1, 0], [1, 1, 1], [0, 1, 0]], + [[0, 0, 0], [0, 1, 0], [0, 0, 0]]], dtype='uint8') + labeled, ncomp = label(img, str_3D) + sizes = ndimage.sum(img, labeled, range(ncomp + 1)) + img_largest_cc = (labeled == np.argmax(sizes)).astype(int) + ## Find coors in each axis for sel_axis in [0, 1, 2]: @@ -15,10 +30,10 @@ def calc_bbox_with_padding(img, perc_pad = 10): other_axes.remove(sel_axis) ## Get img dim in selected axis - dim = img.shape[sel_axis] + dim = img_largest_cc.shape[sel_axis] ## Find bounding box (index of first and last non-zero slices) - nonzero = np.any(img, axis = tuple(other_axes)) + nonzero = np.any(img_largest_cc, axis = tuple(other_axes)) bbox= np.where(nonzero)[0][[0,-1]] ## Add padding From 063a018b59dc950cf2bf4b0ad2969f406b27bf9f Mon Sep 17 00:00:00 2001 From: George Aidins Date: Wed, 13 Sep 2023 16:01:29 -0400 Subject: [PATCH 3/3] Removed duplicate code, extraneous logger line --- niCHARTPipelines/CalculateROIVolume.py | 29 +++++++++++++++++- niCHARTPipelines/MaskImage.py | 42 +------------------------- 2 files changed, 29 insertions(+), 42 deletions(-) diff --git a/niCHARTPipelines/CalculateROIVolume.py b/niCHARTPipelines/CalculateROIVolume.py index 1797993..9a7f944 100644 --- a/niCHARTPipelines/CalculateROIVolume.py +++ b/niCHARTPipelines/CalculateROIVolume.py @@ -24,7 +24,7 @@ def calc_roi_volumes(mrid, in_img_file, label_indices = []): ## Get label indices if label_indices.shape[0] == 0: - logger.warning('Label indices not provided, generating from data') + # logger.warning('Label indices not provided, generating from data') label_indices = u_ind label_names = label_indices.astype(str) @@ -99,3 +99,30 @@ def create_roi_csv(scan_id, in_roi, list_single_roi, map_derived_roi, out_img, o ## Write out csv df_dmuse.to_csv(out_csv, index = False) + +###---------calculate ROI volumes----------- +def extract_roi_masks(in_roi, map_derived_roi, out_pref): + '''Create individual roi masks for single and derived rois + ''' + img_ext_type = '.nii.gz' + + ## Read image + in_nii = nib.load(in_roi) + img_mat = in_nii.get_fdata().astype(int) + + ## Read derived roi map file to a dictionary + roi_dict = {} + with open(map_derived_roi) as roi_map: + reader = csv.reader(roi_map, delimiter=',') + for row in reader: + key = str(row[0]) + val = [int(x) for x in row[2:]] + roi_dict[key] = val + + # Create an individual roi mask for each roi + for i, key in enumerate(roi_dict): + print(i) + key_vals = roi_dict[key] + tmp_mask = np.isin(img_mat, key_vals).astype(int) + out_nii = nib.Nifti1Image(tmp_mask, in_nii.affine, in_nii.header) + nib.save(out_nii, out_pref + '_' + str(key) + img_ext_type) \ No newline at end of file diff --git a/niCHARTPipelines/MaskImage.py b/niCHARTPipelines/MaskImage.py index 6f054a5..161fa2f 100644 --- a/niCHARTPipelines/MaskImage.py +++ b/niCHARTPipelines/MaskImage.py @@ -3,47 +3,7 @@ from scipy import ndimage from scipy.ndimage.measurements import label -## Find bounding box for the foreground values in img, with a given padding percentage -def calc_bbox_with_padding(img, perc_pad = 10): - - img = img.astype('uint8') - - ## Output is the coordinates of the bounding box - bcoors = np.zeros([3,2], dtype=int) - - ## Find the largest connected component - ## INFO: In images with very large FOV DLICV may have small isolated regions in - ## boundaries; so we calculate the bounding box based on the brain, not all - ## foreground voxels - str_3D = np.array([[[0, 0, 0], [0, 1, 0], [0, 0, 0]], - [[0, 1, 0], [1, 1, 1], [0, 1, 0]], - [[0, 0, 0], [0, 1, 0], [0, 0, 0]]], dtype='uint8') - labeled, ncomp = label(img, str_3D) - sizes = ndimage.sum(img, labeled, range(ncomp + 1)) - img_largest_cc = (labeled == np.argmax(sizes)).astype(int) - - ## Find coors in each axis - for sel_axis in [0, 1, 2]: - - ## Get axes other than the selected - other_axes = [0, 1, 2] - other_axes.remove(sel_axis) - - ## Get img dim in selected axis - dim = img_largest_cc.shape[sel_axis] - - ## Find bounding box (index of first and last non-zero slices) - nonzero = np.any(img_largest_cc, axis = tuple(other_axes)) - bbox= np.where(nonzero)[0][[0,-1]] - - ## Add padding - size_pad = int(np.round((bbox[1] - bbox[0]) * perc_pad / 100)) - b_min = int(np.max([0, bbox[0] - size_pad])) - b_max = int(np.min([dim, bbox[1] + size_pad])) - - bcoors[sel_axis, :] = [b_min, b_max] - - return bcoors +from niCHARTPipelines.CombineMasks import calc_bbox_with_padding ###---------mask image-----------