diff --git a/rabies/analysis_pkg/diagnosis_pkg/analysis_QC.py b/rabies/analysis_pkg/diagnosis_pkg/analysis_QC.py index a80db8dd..d204ecd2 100644 --- a/rabies/analysis_pkg/diagnosis_pkg/analysis_QC.py +++ b/rabies/analysis_pkg/diagnosis_pkg/analysis_QC.py @@ -1,8 +1,7 @@ import SimpleITK as sitk import numpy as np import matplotlib.pyplot as plt -import nilearn -from rabies.visualization import otsu_scaling, plot_3d +from rabies.visualization import otsu_scaling, plot_3d, cold_hot from rabies.analysis_pkg.analysis_math import elementwise_spearman, elementwise_corrcoef, dice_coefficient from rabies.utils import recover_3D from rabies.confound_correction_pkg.utils import smooth_image @@ -180,7 +179,7 @@ def masked_plot(fig,axes, img, scaled, mask_img, vmax=None): plot_3d(axes,scaled,fig,vmin=0,vmax=1,cmap='gray', alpha=1, cbar=False, num_slices=6, planes=planes) # resample to match template sitk_img = sitk.Resample(masked, scaled) - cbar_list = plot_3d(axes,sitk_img,fig,vmin=-vmax,vmax=vmax,cmap='cold_hot', alpha=1, cbar=True, threshold=vmax*0.001, num_slices=6, planes=planes) + cbar_list = plot_3d(axes,sitk_img,fig,vmin=-vmax,vmax=vmax,cmap=cold_hot, alpha=1, cbar=True, threshold=vmax*0.001, num_slices=6, planes=planes) return cbar_list diff --git a/rabies/analysis_pkg/diagnosis_pkg/diagnosis_functions.py b/rabies/analysis_pkg/diagnosis_pkg/diagnosis_functions.py index 73ecd24b..fd92b311 100644 --- a/rabies/analysis_pkg/diagnosis_pkg/diagnosis_functions.py +++ b/rabies/analysis_pkg/diagnosis_pkg/diagnosis_functions.py @@ -5,7 +5,6 @@ from rabies.utils import copyInfo_3DImage, recover_3D from rabies.analysis_pkg import analysis_functions import SimpleITK as sitk -import nilearn.plotting from .analysis_QC import masked_plot, threshold_top_percent def compute_spatiotemporal_features(CR_data_dict, sub_maps_data_dict, common_maps_data_dict, analysis_dict, @@ -465,7 +464,7 @@ def scan_diagnosis(CR_data_dict, maps_data_dict, temporal_info, spatial_info, pl fig2, axes2 = plt.subplots(nrows=nrows, ncols=3, figsize=(12*3, 2*nrows)) plt.tight_layout() - from rabies.visualization import otsu_scaling, plot_3d + from rabies.visualization import otsu_scaling, plot_3d, cold_hot axes = axes2[0, :] scaled = otsu_scaling(template_file) @@ -534,7 +533,7 @@ def scan_diagnosis(CR_data_dict, maps_data_dict, temporal_info, spatial_info, pl vector = spatial_info['GS_cov'].flatten() vector.sort() vmax = vector[int(len(vector)*0.95)] - cbar_list = plot_3d(axes, sitk_img, fig2, vmin=-vmax, vmax=vmax, cmap='cold_hot', + cbar_list = plot_3d(axes, sitk_img, fig2, vmin=-vmax, vmax=vmax, cmap=cold_hot, alpha=1, cbar=True, num_slices=6) for cbar in cbar_list: cbar.ax.get_yaxis().labelpad = 20 diff --git a/rabies/confound_correction_pkg/mod_ICA_AROMA/classification_plots.py b/rabies/confound_correction_pkg/mod_ICA_AROMA/classification_plots.py index c653c80b..8385298f 100644 --- a/rabies/confound_correction_pkg/mod_ICA_AROMA/classification_plots.py +++ b/rabies/confound_correction_pkg/mod_ICA_AROMA/classification_plots.py @@ -53,13 +53,13 @@ def classification_plot(myinput, outDir): df3 = pd.DataFrame.from_records([["True", 1., 1., 0., 0.], ["True", 1., 1., 0., 0.], ["True", 1., 1., 0., 0.]]) - df = df.append(df3, ignore_index=True) + df = pd.concat([df, df3], ignore_index=True) tmp = df.loc[df[0] == "False"] if len(tmp) < 3: df3 = pd.DataFrame.from_records([["False", 0., 0., 0., 0.], ["False", 0., 0., 0., 0.], ["False", 0., 0., 0., 0.]]) - df = df.append(df3, ignore_index=True) + df = pd.concat([df, df3], ignore_index=True) # rename columns df = df.rename(index=str, columns={0: 'Motion', diff --git a/rabies/confound_correction_pkg/utils.py b/rabies/confound_correction_pkg/utils.py index 0c31881b..307c3bcb 100644 --- a/rabies/confound_correction_pkg/utils.py +++ b/rabies/confound_correction_pkg/utils.py @@ -642,8 +642,6 @@ def phase_randomized_regressors(regressors_array, frame_mask, TR): def smooth_image(img, fwhm, mask_img): - # apply nilearn's Gaussian smoothing on a SITK image - from nilearn.image.image import _smooth_array from rabies.utils import copyInfo_4DImage, copyInfo_3DImage dim = img.GetDimension() @@ -706,4 +704,94 @@ def sitk_affine_lps(img): affine = np.eye(dim + 1) affine[:dim, :dim] = direction @ spacing affine[:dim, dim] = origin - return affine \ No newline at end of file + return affine + + +""" +Reproduction of Nilearn's nilearn.image.smooth_img smoothing function (version 0.7.1) + +Original reference: +Nilearn developers. Nilearn: Statistical learning for neuroimaging in Python. +URL: https://nilearn.github.io/ +""" +def _smooth_array(arr, affine, fwhm=None, ensure_finite=True, copy=True): + """Smooth images by applying a Gaussian filter. + + Apply a Gaussian filter along the three first dimensions of `arr`. + + Parameters + ---------- + arr : :class:`numpy.ndarray` + 4D array, with image number as last dimension. 3D arrays are also + accepted. + + affine : :class:`numpy.ndarray` + (4, 4) matrix, giving affine transformation for image. (3, 3) matrices + are also accepted (only these coefficients are used). + If `fwhm='fast'`, the affine is not used and can be None. + + fwhm : scalar, :class:`numpy.ndarray`/:obj:`tuple`/:obj:`list`, 'fast' or None, optional + Smoothing strength, as a full-width at half maximum, in millimeters. + If a nonzero scalar is given, width is identical in all 3 directions. + A :class:`numpy.ndarray`, :obj:`tuple`, or :obj:`list` must have 3 elements, + giving the FWHM along each axis. + If any of the elements is zero or None, smoothing is not performed + along that axis. + If `fwhm='fast'`, a fast smoothing will be performed with a filter + [0.2, 1, 0.2] in each direction and a normalisation + to preserve the local average value. + If fwhm is None, no filtering is performed (useful when just removal + of non-finite values is needed). + + ensure_finite : :obj:`bool`, optional + If True, replace every non-finite values (like NaNs) by zero before + filtering. Default=True. + + copy : :obj:`bool`, optional + If True, input array is not modified. True by default: the filtering + is not performed in-place. Default=True. + + Returns + ------- + :class:`numpy.ndarray` + Filtered `arr`. + + Notes + ----- + This function is most efficient with arr in C order. + + """ + # Here, we have to investigate use cases of fwhm. Particularly, if fwhm=0. + # See issue #1537 + if isinstance(fwhm, (int, float)) and (fwhm == 0.0): + import warnings + warnings.warn("The parameter 'fwhm' for smoothing is specified " + "as {0}. Setting it to None " + "(no smoothing will be performed)" + .format(fwhm)) + fwhm = None + if arr.dtype.kind == 'i': + if arr.dtype == np.int64: + arr = arr.astype(np.float64) + else: + arr = arr.astype(np.float32) # We don't need crazy precision. + if copy: + arr = arr.copy() + if ensure_finite: + # SPM tends to put NaNs in the data outside the brain + arr[np.logical_not(np.isfinite(arr))] = 0 + if isinstance(fwhm, str) and (fwhm == 'fast'): + raise NotImplementedError("The 'fast' option for fwhm is not implemented in rabies confound correction. Please specify a numeric value for fwhm or set to None for no smoothing.") + arr = _fast_smooth_array(arr) + elif fwhm is not None: + from scipy.ndimage import gaussian_filter1d + fwhm = np.asarray([fwhm]).ravel() + fwhm = np.asarray([0. if elem is None else elem for elem in fwhm]) + affine = affine[:3, :3] # Keep only the scale part. + fwhm_over_sigma_ratio = np.sqrt(8 * np.log(2)) # FWHM to sigma. + vox_size = np.sqrt(np.sum(affine ** 2, axis=0)) + sigma = fwhm / (fwhm_over_sigma_ratio * vox_size) + for n, s in enumerate(sigma): + if s > 0.0: + gaussian_filter1d(arr, s, output=arr, axis=n) + return arr \ No newline at end of file diff --git a/rabies/preprocess_pkg/preprocess_visual_QC.py b/rabies/preprocess_pkg/preprocess_visual_QC.py index 091171cc..9a43fc1a 100644 --- a/rabies/preprocess_pkg/preprocess_visual_QC.py +++ b/rabies/preprocess_pkg/preprocess_visual_QC.py @@ -65,7 +65,6 @@ def template_info(anat_template, opts, out_dir,figure_format): import SimpleITK as sitk # set default threader to platform to avoid freezing with MultiProc https://github.com/SimpleITK/SimpleITK/issues/1239 sitk.ProcessObject_SetGlobalDefaultThreader('Platform') - from nilearn import plotting import matplotlib.pyplot as plt from rabies.visualization import plot_3d, otsu_scaling os.makedirs(out_dir, exist_ok=True) @@ -99,7 +98,6 @@ def template_masking(template, mask, out_dir, figure_format): import SimpleITK as sitk # set default threader to platform to avoid freezing with MultiProc https://github.com/SimpleITK/SimpleITK/issues/1239 sitk.ProcessObject_SetGlobalDefaultThreader('Platform') - from nilearn import plotting import matplotlib.pyplot as plt from rabies.visualization import plot_3d, otsu_scaling diff --git a/rabies/visualization.py b/rabies/visualization.py index 42ede81f..d06a8b30 100644 --- a/rabies/visualization.py +++ b/rabies/visualization.py @@ -18,6 +18,46 @@ "savefig.facecolor": "black", "savefig.edgecolor": "black"}) +""" +Standalone reproduction of Nilearn's 'cold_hot' colormap. + +Original reference: +Nilearn developers. Nilearn: Statistical learning for neuroimaging in Python. +URL: https://nilearn.github.io/ +""" +from matplotlib.colors import LinearSegmentedColormap +# Segment data extracted from the actual Nilearn 'cold_hot' +_cold_hot_segments = { + 'red': [ + (0.0, 1.0, 1.0), + (0.126984, 0.0, 0.0), + (0.5, 0.0, 0.0), + (0.5, 0.0416, 0.0416), + (0.682540, 1.0, 1.0), + (1.0, 1.0, 1.0), + ], + 'green': [ + (0.0, 1.0, 1.0), + (0.126984, 1.0, 1.0), + (0.317461, 0.0, 0.0), + (0.5, 0.0, 0.0), + (0.5, 0.0, 0.0), + (0.682540, 0.0, 0.0), + (0.873016, 1.0, 1.0), + (1.0, 1.0, 1.0), + ], + 'blue': [ + (0.0, 1.0, 1.0), + (0.317461, 1.0, 1.0), + (0.5, 0.0416, 0.0416), + (0.5, 0.0, 0.0), + (0.873016, 0.0, 0.0), + (1.0, 1.0, 1.0), + ] +} +# Create the LinearSegmentedColormap +cold_hot = LinearSegmentedColormap('cold_hot', _cold_hot_segments, N=256) + def otsu_scaling(input_image): from skimage.filters import threshold_multiotsu diff --git a/rabies_environment.dev.yml b/rabies_environment.dev.yml index ecf40d61..98496049 100644 --- a/rabies_environment.dev.yml +++ b/rabies_environment.dev.yml @@ -9,7 +9,6 @@ dependencies: - future - matplotlib - nibabel - - nilearn - nipype - numpy - pandas diff --git a/rabies_environment.yml b/rabies_environment.yml index a5bc2715..6889731b 100644 --- a/rabies_environment.yml +++ b/rabies_environment.yml @@ -9,7 +9,6 @@ dependencies: - future - matplotlib=3.3.4 - nibabel=3.2.1 - - nilearn=0.7.1 - nipype=1.10.0 - numpy=1.26.4 - pandas=1.2.4 diff --git a/setup.py b/setup.py index 72a71e78..d687aefb 100644 --- a/setup.py +++ b/setup.py @@ -17,26 +17,24 @@ URL = 'https://github.com/CoBrALab/RABIES' EMAIL = 'contact@cobralab.ca' AUTHOR = 'CoBrALab' -REQUIRES_PYTHON = '>=3.9.0' +REQUIRES_PYTHON = '>=3.9.0, <3.13' # the package was tested up to 3.12 -# What packages are required for this module to be executed? REQUIRED = [ # 'requests', 'maya', 'records', 'pip>=23.0', - 'matplotlib>=3.3.4', - 'nibabel>=3.2.1', - 'nilearn>=0.7.1', - 'nipype>=1.6.1', - 'numpy>=1.20.1', - 'pandas>=1.2.4', - 'pathos>=0.2.7', - 'pybids==0.16.3', - 'scikit-learn>=0.24.1', - 'scikit-image>=0.18.2', - 'scipy>=1.8.1', - 'seaborn>=0.11.1', + 'matplotlib>=3.3.4', # tested up to 3.10.8 + 'nibabel>=3.2.1', # tested up to 5.4.0 + 'nipype>=1.6.1', # tested up to 1.11.0 + 'numpy>=1.20.1', # tested up to 2.4.3 + 'pandas>=1.2.4', # tested up to 3.0.1 + 'pathos>=0.2.7', # tested up to 0.3.5 + 'pybids==0.16.3', # only tested with 0.16.3, newer versions have some issues with the BIDS layout of our data + 'scikit-learn>=0.24.1', # tested up to 1.8.0 + 'scikit-image>=0.18.2', # tested up to 0.24.0 + 'scipy>=1.8.1', # tested up to 1.17.1 + 'seaborn>=0.11.1', # tested up to 0.13.2 'simpleitk>=2.5.0', - 'qbatch==2.3', + 'qbatch==2.3', # tested only with 2.3 'networkx<3', 'traits<7.0', 'tqdm',