diff --git a/environment_cpu.yml b/archive/environment_cpu.yml similarity index 100% rename from environment_cpu.yml rename to archive/environment_cpu.yml diff --git a/environment_gpu.yml b/archive/environment_gpu.yml similarity index 100% rename from environment_gpu.yml rename to archive/environment_gpu.yml diff --git a/docs/README.md b/docs/README.md index a02c488..1d33f6f 100644 --- a/docs/README.md +++ b/docs/README.md @@ -8,35 +8,35 @@ Brigham and Women's Hospital (Harvard Medical School). Table of Contents ================= - * [Table of Contents](#table-of-contents) - * [Segmenting diffusion MRI brain](#segmenting-diffusion-mri-brain) - * [Citation](#citation) - * [Dependencies](#dependencies) - * [Installation](#installation) - * [1. Python 3](#1-python-3) - * [2. Conda environment](#2-conda-environment) - * [CPU only](#cpu-only) - * [GPU support](#gpu-support) - * [3. CUDA environment](#3-cuda-environment) - * [4. Download models](#4-download-models) - * [5. mrtrix](#5-mrtrix) - * [Singularity container](#singularity-container) - * [Running the pipeline](#running-the-pipeline) - * [Prediction](#prediction) - * [Training](#training) - * [1. B0 extraction](#1-b0-extraction) - * [2. Registration](#2-registration) - * [3. Preprocessing](#3-preprocessing) - * [4. Deep learning](#4-deep-learning) - * [Method](#method) - * [1. Model Architecture](#1-model-architecture) - * [2. Multi View Aggregation:](#2-multi-view-aggregation) - * [3. Clean up](#3-clean-up) - * [Troubleshooting](#troubleshooting) - * [Issues](#issues) - * [Reference](#reference) - -Table of contents created by [gh-md-toc](https://github.com/ekalinin/github-markdown-toc) +* [Table of Contents](#table-of-contents) +* [Segmenting diffusion MRI brain](#segmenting-diffusion-mri-brain) +* [Citation](#citation) +* [Dependencies](#dependencies) +* [Installation](#installation) + * [1. Python 3](#1-python-3) + * [2. Conda environment](#2-conda-environment) + * [3. ANTs](#3-ants) + * [4. CUDA environment](#4-cuda-environment) + * [5. Download models](#5-download-models) + * [6. mrtrix](#6-mrtrix) + * [7. FSL](#7-fsl) +* [Singularity container](#singularity-container) +* [Running the pipeline](#running-the-pipeline) + * [Prediction](#prediction) + * [Training](#training) + * [1. B0 extraction](#1-b0-extraction) + * [2. Registration](#2-registration) + * [3. Preprocessing](#3-preprocessing) + * [4. Deep learning](#4-deep-learning) +* [Method](#method) + * [1. Model Architecture](#1-model-architecture) + * [2. Multi View Aggregation:](#2-multi-view-aggregation) + * [3. Clean up](#3-clean-up) +* [Troubleshooting](#troubleshooting) +* [Issues](#issues) +* [Reference](#reference) + + # Segmenting diffusion MRI brain @@ -71,7 +71,7 @@ NeuroImage 186 (2019): 713-727. # Dependencies -* python 3 +* Python 3 * ANTs # Installation @@ -86,47 +86,71 @@ Activate the conda environment: source ~/miniconda3/bin/activate # should introduce '(base)' in front of each line -The software is written intelligently to work with or without GPU support. Two *yml* environment files -are provided to facilitate creation of `dmri_seg` conda environment. +The software is written intelligently to work with or without GPU support. +
Old conda environment + +
+ +We had provided two [*yml* files](../archive/) to facilitate the creation of `dmri_seg` conda environment. + +* In the initial stage of this program development in 2020, the [CPU environment](../archive/environment_cpu.yml) built and worked fine. + But it does not work anymore. +* The [GPU environment](../archive/environment_gpu.yml) works fine on older GPUs. + But in 2023, it did not produce a valid mask on modern GPUs e.g. NVIDIA RTX 4000, A6000, A100. + +
+ +In 2024, we recommend the below step-by-step instructions to build an environment that + successfully generated valid masks on both GPU and CPU devices. ## 2. Conda environment - -### CPU only - conda env create -f environment_cpu.yml +Step-by-step instructions: +``` +conda create -y -n dmri_seg python=3.9 +conda activate dmri_seg +pip install tensorflow==2.11 +conda install -y anaconda::cudnn conda-forge::gputil +pip install nibabel scikit-image git+https://github.com/pnlbwh/conversion.git +``` -### GPU support +## 3. ANTs - conda env create -f environment_gpu.yml - -Finally, activate the conda environment using: +You need to have a working `antsRegistrationSyNQuick.sh` in your `PATH`. +You can build ANTs following their [GitHub](https://github.com/ANTsX/ANTs). +Or you can use our pre-built ANTs: - conda activate dmri_seg - +> conda install -y pnlbwh::ants -If you run into any error, please see [Troubleshooting](#troubleshooting). +Either way, you need to set the environment variable `ANTSPATH`. For the latter method, it is: +> export ANTSPATH=${CONDA_PREFIX}/bin -## 3. CUDA environment +## 4. CUDA environment -When you have GPU support, provided that you used `environment_gpu.yml` for creating conda environment, -you should set environment variables in order to run and write CUDA enabled programs. -The NVIDIA graphics driver and CUDA compilier are already installed on machines that support CUDA. +To run this program on GPU, you must set `LD_LIBRARY_PATH`: -If you use bash, add the following lines to the bottom of your `~/.bashrc` file: +> export LD_LIBRARY_PATH=${CONDA_PREFIX}/lib/:${LD_LIBRARY_PATH} - # add cuda tools to command path - export PATH=/usr/local/cuda-9.1/bin:${PATH} +Your NVIDIA driver should be compatible with CUDA. - # add the CUDA binary and library directory to your LD_LIBRARY_PATH - export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/cuda-9.1/lib64 +
System's CUDA +
-Open a new terminal for the changes to take effect. +In some servers, a matched CUDA may be available in `/usr/local/cuda-*/` +directory. You can also use that CUDA instead of `${CONDA_PREFIX}/lib/`: + + # add cuda tools to your PATH + export PATH=/usr/local/cuda-9.1/bin:${PATH} + + # add cuda libraries to your LD_LIBRARY_PATH + export LD_LIBRARY_PATH=/usr/local/cuda-9.1/lib64:${LD_LIBRARY_PATH} +
-## 4. Download models +## 5. Download models Download model architecture, weights and IIT mean b0 template from https://github.com/pnlbwh/CNN-Diffusion-MRIBrain-Segmentation/releases as follows: @@ -141,7 +165,7 @@ https://github.com/pnlbwh/CNN-Diffusion-MRIBrain-Segmentation/releases as follow They will be extracted to `CNN-Diffusion-MRIBrain-Segmentation/model_folder` directory. -## 5. mrtrix +## 6. mrtrix This software makes use of only one binary from mrtrix, `maskfilter`. If you already have mrtrix installed on your system, you can use it. Moreover, you can follow [this](https://mrtrix.readthedocs.io/en/latest/installation/linux_install.html) instruction to install mrtrix. @@ -151,6 +175,19 @@ use a Python translated version of the above binary. See [Clean up](#3-clean-up) for details. +## 7. FSL + +This software marginally depends on FSL in the sense that you need FSL's `slicesdir` executable +to generate snapshots of image-mask for being able to QC it. Hence, we do not recommend +sourcing the entire FSL environment. Instead, you should just put `slicesdir` in your PATH: + +> export PATH=/path/to/fsl-6.0.7/share/fsl/bin:$PATH + +In additon, you should set: + +> export FSLOUTPUTTYPE=NIFTI_GZ + + # Singularity container The dependencies in [environment_cpu.yml](../environment_cpu.yml) have been found to be volatile. So we have archived @@ -185,7 +222,9 @@ Prediction refers to creation of masks based on pre-trained model. This is the c pipeline/dwi_masking.py -i dwi_list.txt -f model_folder pipeline/dwi_masking.py -i dwi_list.txt -f model_folder -nproc 16 pipeline/dwi_masking.py -i b0_list.txt -f model_folder - + pipeline/dwi_masking.py -i b0_list.txt -f model_folder -qc 1 + + * `dwi_list.txt` and `b0_list.txt` must contain full paths to diffusion and b0 volumes respectively * Each created mask is saved in the directory of input volume with name `dwib0_{PREFIX}-multi_BrainMask.nii.gz` diff --git a/pipeline/dwi_masking.py b/pipeline/dwi_masking.py index 4f168e5..e8fd11e 100755 --- a/pipeline/dwi_masking.py +++ b/pipeline/dwi_masking.py @@ -1,12 +1,6 @@ #!/usr/bin/env python from __future__ import division -# ----------------------------------------------------------------- -# Author: Senthil Palanivelu, Tashrif Billah -# Written: 01/22/2020 -# Last Updated: 02/28/2020 -# Purpose: Pipeline for diffusion brain masking -# ----------------------------------------------------------------- """ pipeline.py @@ -15,21 +9,30 @@ 02) Checks if the Image axis is in the correct order for *.nhdr and *.nrrd file 03) Extracts b0 Image 04) Converts nhdr to nii.gz -05) Applys Rigid-Body tranformation to standard MNI space using +05) Applies rigid-body tranformation to standard MNI space using 06) Normalize the Image by 99th percentile -07) Neural network brain mask prediction across the 3 principal axis -08) Performs Multi View Aggregation -10) Applys Inverse tranformation -10) Cleaning +07) Predicts neural network brain mask across the 3 principal axes +08) Performs multi-view aggregation +10) Applies inverse tranformation +10) Cleans holes """ - # pylint: disable=invalid-name import os from os import path import webbrowser -os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2' # Suppress tensor flow message +import multiprocessing as mp +import sys +from glob import glob +import subprocess +import argparse +import datetime +import pathlib +import nibabel as nib +import numpy as np +# Suppress tensor flow message +os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2' # Set CUDA_DEVICE_ORDER so the IDs assigned by CUDA match those from nvidia-smi os.environ["CUDA_DEVICE_ORDER"] = "PCI_BUS_ID" @@ -37,27 +40,26 @@ # Get the first available GPU try: import GPUtil + DEVICE_ID_LIST = [g.id for g in GPUtil.getGPUs()] DEVICE_ID = DEVICE_ID_LIST[0] - - - CUDA_VISIBLE_DEVICES=os.getenv('CUDA_VISIBLE_DEVICES') + + CUDA_VISIBLE_DEVICES = os.getenv('CUDA_VISIBLE_DEVICES') if CUDA_VISIBLE_DEVICES: # prioritize external definition if int(CUDA_VISIBLE_DEVICES) in DEVICE_ID_LIST: pass else: # define it internally - CUDA_VISIBLE_DEVICES=DEVICE_ID + CUDA_VISIBLE_DEVICES = DEVICE_ID else: # define it internally - CUDA_VISIBLE_DEVICES=DEVICE_ID - + CUDA_VISIBLE_DEVICES = DEVICE_ID + os.environ["CUDA_VISIBLE_DEVICES"] = str(CUDA_VISIBLE_DEVICES) # setting of CUDA_VISIBLE_DEVICES also masks out all other GPUs - - print ("Use GPU device #", CUDA_VISIBLE_DEVICES) + print("Use GPU device #", CUDA_VISIBLE_DEVICES) except: os.environ["CUDA_VISIBLE_DEVICES"] = "-1" @@ -65,36 +67,39 @@ import warnings + with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=FutureWarning) import tensorflow as tf - -import multiprocessing as mp -import re -import sys -from glob import glob -import subprocess -import argparse, textwrap -import datetime -import pathlib -import nibabel as nib -import numpy as np -from keras.models import load_model -from keras.models import model_from_json -from multiprocessing import Process, Manager, Value, Pool -from time import sleep -import keras -from keras import losses -from keras.models import Model -from keras.layers import Input, merge, concatenate, Conv2D, MaxPooling2D, \ - Activation, UpSampling2D, Dropout, Conv2DTranspose, add, multiply -from keras.layers.normalization import BatchNormalization as bn -from keras.callbacks import ModelCheckpoint, TensorBoard -from keras.optimizers import RMSprop -from keras import regularizers -from keras import backend as K -from keras.optimizers import Adam -from keras.callbacks import ModelCheckpoint + try: + from keras.models import model_from_json + from keras import backend as K + from keras.optimizers import Adam + except ImportError: + from tensorflow.keras.models import model_from_json + from tensorflow.keras import backend as K + from tensorflow.keras.optimizers import Adam + + # check version of tf and if 1.12 or less use tf.logging.set_verbosity(tf.logging.ERROR) + if int(tf.__version__.split('.')[0]) <= 1 and int(tf.__version__.split('.')[1]) <= 12: + tf.logging.set_verbosity(tf.logging.ERROR) + # Configure for dynamic GPU memory usage + config = tf.ConfigProto() + config.gpu_options.allow_growth = True + config.log_device_placement = False + sess = tf.Session(config=config) + K.set_session(sess) + else: + tf.compat.v1.logging.set_verbosity(tf.compat.v1.logging.ERROR) + gpus = tf.config.experimental.list_physical_devices('GPU') + if gpus: + try: + # Currently, memory growth needs to be the same across GPUs + for gpu in gpus: + tf.config.experimental.set_memory_growth(gpu, True) + except RuntimeError as e: + # Memory growth must be set before GPUs have been initialized + print(e) # suffixes SUFFIX_NIFTI = "nii" @@ -122,7 +127,7 @@ def predict_mask(input_file, trained_folder, view='default'): returns the neural network predicted filename which is stored in disk in 3d numpy array *.npy format """ - print ("Loading " + view + " model from disk...") + print("Loading " + view + " model from disk...") smooth = 1. def dice_coef(y_true, y_pred): @@ -144,16 +149,22 @@ def neg_dice_coef_loss(y_true, y_pred): loaded_model_json = json_file.read() json_file.close() loaded_model = model_from_json(loaded_model_json) - + # load weights into new model - optimal_model= glob(trained_folder + '/weights-' + view + '-improvement-*.h5')[-1] + optimal_model = glob(trained_folder + '/weights-' + view + '-improvement-*.h5')[-1] loaded_model.load_weights(optimal_model) - # evaluate loaded model on test data - loaded_model.compile(optimizer=Adam(lr=1e-5), - loss={'final_op': dice_coef_loss, - 'xfinal_op': neg_dice_coef_loss, - 'res_1_final_op': 'mse'}) + # check if tf2 or tf1, if tf 1 use lr instead of learning_rate + if int(tf.__version__.split('.')[0]) <= 1: + loaded_model.compile(optimizer=Adam(lr=1e-5), + loss={'final_op': dice_coef_loss, + 'xfinal_op': neg_dice_coef_loss, + 'res_1_final_op': 'mse'}) + else: + loaded_model.compile(optimizer=Adam(learning_rate=1e-5), + loss={'final_op': dice_coef_loss, + 'xfinal_op': neg_dice_coef_loss, + 'res_1_final_op': 'mse'}) case_name = path.basename(input_file) output_name = case_name[:len(case_name) - (len(SUFFIX_NIFTI_GZ) + 1)] + '-' + view + '-mask.npy' @@ -169,7 +180,6 @@ def neg_dice_coef_loss(y_true, y_pred): def multi_view_fast(sagittal_SO, coronal_SO, axial_SO, input_file): - x = np.load(sagittal_SO) y = np.load(coronal_SO) z = np.load(axial_SO) @@ -202,7 +212,7 @@ def multi_view_fast(sagittal_SO, coronal_SO, axial_SO, input_file): return output_file -def normalize(b0_resampled, percentile, data_n): +def normalize(b0_resampled, percentile): """ Intensity based segmentation of MR images is hampered by radio frerquency field inhomogeneity causing intensity variation. The intensity range is typically @@ -220,25 +230,24 @@ def normalize(b0_resampled, percentile, data_n): output_file : str Normalized by 99th percentile filename which is stored in disk """ - print ("Normalizing input data") + print("Normalizing input data") input_file = b0_resampled case_name = path.basename(input_file) output_name = case_name[:len(case_name) - (len(SUFFIX_NIFTI_GZ) + 1)] + '-normalized.nii.gz' output_file = path.join(path.dirname(input_file), output_name) img = nib.load(b0_resampled) - imgU16 = img.get_data().astype(np.float32) + imgU16 = img.get_fdata().astype(np.float32) p = np.percentile(imgU16, percentile) data = imgU16 / p data[data > 1] = 1 data[data < 0] = 0 image_dwi = nib.Nifti1Image(data, img.affine, img.header) nib.save(image_dwi, output_file) - data_n.append(output_file) + return output_file def save_nifti(fname, data, affine=None, hdr=None): - hdr.set_data_dtype('int16') result_img = nib.Nifti1Image(data, affine, header=hdr) result_img.to_filename(fname) @@ -311,24 +320,23 @@ def npy_to_nifti(b0_normalized_cases, cases_mask_arr, sub_name, view='default', if args.filter: print('Cleaning up ', CNN_output_file) - - if args.filter=='mrtrix': + + if args.filter == 'mrtrix': mask_filter = "maskfilter -force " + CNN_output_file + " -scale 2 clean " + output_mask_filtered - elif args.filter=='scipy': - mask_filter = path.join(path.dirname(__file__),'../src/maskfilter.py') + f' {CNN_output_file} 2 {output_mask_filtered}' - + elif args.filter == 'scipy': + mask_filter = path.join(path.dirname(__file__), '../src/maskfilter.py') + f' {CNN_output_file} 2 {output_mask_filtered}' + process = subprocess.Popen(mask_filter.split(), stdout=subprocess.PIPE) output, error = process.communicate() else: - output_mask_filtered= CNN_output_file - + output_mask_filtered = CNN_output_file print(output_mask_filtered) img = nib.load(output_mask_filtered) data_dwi = nib.load(sub_name[i]) - imgU16 = img.get_data().astype(np.uint8) + imgU16 = img.get_fdata().astype(np.uint8) brain_mask_file = subject_name[:len(subject_name) - (len(format) + 1)] + '-' + view + '_BrainMask.nii.gz' brain_mask_final = path.join(output_dir, brain_mask_file) @@ -340,7 +348,7 @@ def npy_to_nifti(b0_normalized_cases, cases_mask_arr, sub_name, view='default', def clear(directory): - print ("Cleaning files ...") + print("Cleaning files ...") bin_a = 'cases_' + str(os.getpid()) + '_binary_a' bin_s = 'cases_' + str(os.getpid()) + '_binary_s' @@ -352,11 +360,11 @@ def clear(directory): filename.endswith('-thresholded.nii.gz') | filename.endswith('-inverse.mat') | \ filename.endswith('-Warped.nii.gz') | filename.endswith('-0GenericAffine.mat') | \ filename.endswith('_affinedMask.nii.gz') | filename.endswith('_originalMask.nii.gz') | \ - filename.endswith('multi-mask.nii.gz') | filename.endswith('-mask-inverse.nii.gz') | \ + filename.endswith('multi-mask.nii.gz') | filename.endswith('-mask-inverse.nii.gz') | \ filename.endswith('-InverseWarped.nii.gz') | filename.endswith('-FilteredMask.nii.gz') | \ filename.endswith(bin_a) | filename.endswith(bin_c) | filename.endswith(bin_s) | \ filename.endswith('_FilteredMask.nii.gz') | filename.endswith('-normalized.nii.gz') | filename.endswith('-filled.nii.gz'): - os.unlink(directory + '/' + filename) + os.unlink(directory + '/' + filename) def split(cases_file, case_arr, view='default'): @@ -375,7 +383,6 @@ def split(cases_file, case_arr, view='default'): Contains the predicted mask filename of all the cases which is stored in disk in *.npy format """ - count = 0 start = 0 end = start + 256 @@ -390,7 +397,7 @@ def split(cases_file, case_arr, view='default'): elif view == 'axial': casex = np.swapaxes(casex, 0, 2) input_file = str(case_arr[i]) - output_file = input_file[:len(input_file) - (len(SUFFIX_NHDR) + 1)] + '-' + view +'_SO.npy' + output_file = input_file[:len(input_file) - (len(SUFFIX_NHDR) + 1)] + '-' + view + '_SO.npy' predict_mask.append(output_file) np.save(output_file, casex) start = end @@ -399,8 +406,7 @@ def split(cases_file, case_arr, view='default'): return predict_mask -def ANTS_rigid_body_trans(b0_nii, result, reference=None): - +def ANTS_rigid_body_trans(b0_nii, reference=None): print("Performing ants rigid body transformation...") input_file = b0_nii case_name = path.basename(input_file) @@ -416,15 +422,11 @@ def ANTS_rigid_body_trans(b0_nii, result, reference=None): output_name = case_name[:len(case_name) - (len(SUFFIX_NIFTI_GZ) + 1)] + '-Warped.nii.gz' transformed_file = path.join(path.dirname(input_file), output_name) - result.append((transformed_file, omat_file)) + return (transformed_file, omat_file) def ANTS_inverse_transform(predicted_mask, reference, omat='default'): - #print "Mask file = ", predicted_mask - #print "Reference = ", reference - #print "omat = ", omat - print("Performing ants inverse transform...") input_file = predicted_mask case_name = path.basename(input_file) @@ -433,7 +435,7 @@ def ANTS_inverse_transform(predicted_mask, reference, omat='default'): # reference is the original b0 volume apply_inverse_trans = "antsApplyTransforms -d 3 -i " + predicted_mask + " -r " + reference + " -o " \ - + output_file + " --transform [" + omat + ",1]" + + output_file + " --transform [" + omat + ",1]" output2 = subprocess.check_output(apply_inverse_trans, shell=True) return output_file @@ -441,7 +443,7 @@ def ANTS_inverse_transform(predicted_mask, reference, omat='default'): def str2bool(v): if isinstance(v, bool): - return v + return v if v.lower() in ('yes', 'true', 't', 'y', '1'): return True elif v.lower() in ('no', 'false', 'f', 'n', '0'): @@ -451,43 +453,40 @@ def str2bool(v): def list_masks(mask_list, view='default'): - for i in range(0, len(mask_list)): - print (view + " Mask file = ", mask_list[i]) + print(view + " Mask file = ", mask_list[i]) -def pre_process(input_file, target_list, b0_threshold=50.): - +def pre_process(input_file, b0_threshold=50.): from conversion import nifti_write, read_bvals - from subprocess import Popen if path.isfile(input_file): # convert NRRD/NHDR to NIFIT as the first step # extract bse.py from just NIFTI later if input_file.endswith(SUFFIX_NRRD) | input_file.endswith(SUFFIX_NHDR): - inPrefix= input_file.split('.')[0] + inPrefix = input_file.split('.')[0] nifti_write(input_file) - input_file= inPrefix+ '.nii.gz' + input_file = inPrefix + '.nii.gz' + + inPrefix = input_file.split('.nii')[0] + b0_nii = path.join(inPrefix + '_bse.nii.gz') - inPrefix= input_file.split('.nii')[0] - b0_nii= path.join(inPrefix+ '_bse.nii.gz') - - dwi= nib.load(input_file) + dwi = nib.load(input_file) - if len(dwi.shape)>3: + if len(dwi.shape) > 3: print("Extracting b0 volume...") - bvals= np.array(read_bvals(input_file.split('.nii')[0]+ '.bval')) - where_b0= np.where(bvals <= b0_threshold)[0] - b0= dwi.get_data()[...,where_b0].mean(-1) + bvals = np.array(read_bvals(input_file.split('.nii')[0] + '.bval')) + where_b0 = np.where(bvals <= b0_threshold)[0] + b0 = dwi.get_fdata()[..., where_b0].mean(-1) else: print("Loading b0 volume...") - b0= dwi.get_fdata() + b0 = dwi.get_fdata() - np.nan_to_num(b0).clip(min= 0., out= b0) - nib.Nifti1Image(b0, affine= dwi.affine, header= dwi.header).to_filename(b0_nii) + np.nan_to_num(b0).clip(min=0., out=b0) + nib.Nifti1Image(b0, affine=dwi.affine, header=dwi.header).to_filename(b0_nii) - target_list.append((b0_nii)) + return b0_nii else: print("File not found ", input_file) @@ -523,24 +522,24 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): dir_bak = os.getcwd() os.chdir(tmp_path) - process= subprocess.Popen(final, shell=True) + process = subprocess.Popen(final, shell=True) process.wait() os.chdir(dir_bak) mask_folder = os.path.join(tmp_path, 'slicesdir') mask_newfolder = os.path.join(tmp_path, 'slicesdir_' + view) if os.path.exists(mask_newfolder): - process = subprocess.Popen('rm -rf '+ mask_newfolder, shell=True) + process = subprocess.Popen('rm -rf ' + mask_newfolder, shell=True) process.wait() process = subprocess.Popen('mv ' + mask_folder + " " + mask_newfolder, shell=True) process.wait() - + if __name__ == '__main__': start_total_time = datetime.datetime.now() - # parser module for input arguments + parser = argparse.ArgumentParser() parser.add_argument('-i', action='store', dest='dwi', type=str, @@ -563,14 +562,14 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): parser.add_argument("-qc", type=str2bool, dest='snap', nargs='?', const=True, default=False, - help="open snapshots in your web browser (yes/true/y/1)") + help="advanced option to take snapshots and open them in your web browser (yes/true/y/1)") - parser.add_argument('-p', type=int, dest='percentile', default=99, help='''The percentile of image + parser.add_argument('-p', type=int, dest='percentile', default=99, help='''percentile of image intensity value to be used as a threshold for normalizing a b0 image to [0,1]''') - parser.add_argument('-nproc', type=int, dest='cr', default=8, help='number of processes to use') - - parser.add_argument('-filter', choices=['scipy','mrtrix'], help='''perform morphological operation on the + parser.add_argument('-nproc', type=int, dest='nproc', default=1, help='number of processes to use') + + parser.add_argument('-filter', choices=['scipy', 'mrtrix'], help='''perform morphological operation on the CNN generated mask to clean up holes and islands, can be done through a provided script (scipy) or MRtrix3 maskfilter (mrtrix)''') @@ -580,17 +579,16 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): parser.print_help() parser.error('too few arguments') sys.exit(0) - except SystemExit: sys.exit(0) if args.dwi: f = pathlib.Path(args.dwi) if f.exists(): - print ("File exist") + print("File exist") filename = args.dwi else: - print ("File not found") + print("File not found") sys.exit(1) # Input caselist.txt @@ -598,11 +596,8 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): with open(filename) as f: case_arr = f.read().splitlines() - TXT_file = path.basename(filename) - #print(TXT_file) - unique = TXT_file[:len(TXT_file) - (len(SUFFIX_TXT)+1)] - #print(unique) + unique = TXT_file[:len(TXT_file) - (len(SUFFIX_TXT) + 1)] storage = path.dirname(case_arr[0]) tmp_path = storage + '/' trained_model_folder = args.model_folder.rstrip('/') @@ -620,83 +615,76 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): y_dim = 256 z_dim = 256 transformed_cases = [] - """ - Enable Multi core Processing for pre processing - manager provide a way to create data which can be shared between different processes - """ - with Manager() as manager: - target_list = manager.list() - omat_list = [] - jobs = [] - for i in range(0,len(case_arr)): - p_process = mp.Process(target=pre_process, args=(case_arr[i], - target_list)) - p_process.start() - jobs.append(p_process) - - for process in jobs: - process.join() - - target_list = list(target_list) - """ - Enable Multi core Processing for ANTS Registration - manager provide a way to create data which can be shared between different processes - """ - with Manager() as manager: - result = manager.list() - ants_jobs = [] - for i in range(0, len(target_list)): - p_ants = mp.Process(target=ANTS_rigid_body_trans, args=(target_list[i], - result, reference)) - ants_jobs.append(p_ants) - p_ants.start() - - for process in ants_jobs: - process.join() - - result = list(result) - - for subject_ANTS in result: - transformed_cases.append(subject_ANTS[0]) - omat_list.append(subject_ANTS[1]) - - with Manager() as manager: - data_n = manager.list() - norm_jobs = [] - for i in range(0, len(target_list)): - p_norm = mp.Process(target=normalize, args=(transformed_cases[i], - args.percentile, data_n)) - norm_jobs.append(p_norm) - p_norm.start() - - for process in norm_jobs: - process.join() - - data_n = list(data_n) - - + + if args.nproc==1: + + target_list=[] + for case in case_arr: + target_list.append(pre_process(case)) + + result=[] + for target in target_list: + result.append(ANTS_rigid_body_trans(target,reference)) + + data_n=[] + for transformed_case, _ in result: + data_n.append(normalize(transformed_case, args.percentile)) + + else: + with mp.Pool(processes=args.nproc) as pool: + res=[] + for case in case_arr: + res.append(pool.apply_async(pre_process, (case,))) + + target_list=[r.get() for r in res] + pool.close() + pool.join() + + + with mp.Pool(processes=args.nproc) as pool: + res=[] + for target in target_list: + res.append(pool.apply_async(ANTS_rigid_body_trans, (target, reference,))) + + result=[r.get() for r in res] + pool.close() + pool.join() + + + with mp.Pool(processes=args.nproc) as pool: + res=[] + for transformed_case, _ in result: + res.append(pool.apply_async(normalize, (transformed_case, args.percentile,))) + + data_n=[r.get() for r in res] + pool.close() + pool.join() + + count = 0 for b0_nifti in data_n: img = nib.load(b0_nifti) - imgU16_sagittal = img.get_data().astype(np.float32) # sagittal view - imgU16_coronal = np.swapaxes(imgU16_sagittal, 0, 1) # coronal view - imgU16_axial = np.swapaxes(imgU16_sagittal, 0, 2) # Axial view + # sagittal view + imgU16_sagittal = img.get_fdata().astype(np.float32) + # coronal view + imgU16_coronal = np.swapaxes(imgU16_sagittal, 0, 1) + # axial view + imgU16_axial = np.swapaxes(imgU16_sagittal, 0, 2) imgU16_sagittal.tofile(f_handle_s) imgU16_coronal.tofile(f_handle_c) imgU16_axial.tofile(f_handle_a) - - print ("Case completed = ", count) count += 1 + print("Case completed = ", count) f_handle_s.close() f_handle_c.close() f_handle_a.close() - print ("Merging npy files...") - cases_file_s = storage + '/'+ unique + '_' + str(os.getpid()) + '-casefile-sagittal.npy' - cases_file_c = storage + '/'+ unique + '_' + str(os.getpid()) + '-casefile-coronal.npy' - cases_file_a = storage + '/'+ unique + '_' + str(os.getpid()) + '-casefile-axial.npy' + print("Merging npy files...") + cases_file_s = storage + '/' + unique + '_' + str(os.getpid()) + '-casefile-sagittal.npy' + cases_file_c = storage + '/' + unique + '_' + str(os.getpid()) + '-casefile-coronal.npy' + cases_file_a = storage + '/' + unique + '_' + str(os.getpid()) + '-casefile-axial.npy' merged_dwi_list = [] merged_dwi_list.append(cases_file_s) @@ -707,7 +695,7 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): merge_c = np.memmap(binary_file_c, dtype=np.float32, mode='r+', shape=(256 * len(target_list), y_dim, z_dim)) merge_a = np.memmap(binary_file_a, dtype=np.float32, mode='r+', shape=(256 * len(target_list), y_dim, z_dim)) - print ("Saving data to disk...") + print("Saving data to disk...") np.save(cases_file_s, merge_s) np.save(cases_file_c, merge_c) np.save(cases_file_a, merge_a) @@ -725,14 +713,14 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): remove_string(registered_file, target_file, "-Warped") with open(target_file) as f: - newText=f.read().replace('.nii.gz', '-0GenericAffine.mat') + newText = f.read().replace('.nii.gz', '-0GenericAffine.mat') with open(mat_file, "w") as f: f.write(newText) end_preprocessing_time = datetime.datetime.now() total_preprocessing_time = end_preprocessing_time - start_total_time - print ("Pre-Processing Time Taken : ", round(int(total_preprocessing_time.seconds)/60, 2), " min") + print("Pre-Processing Time Taken : ", round(int(total_preprocessing_time.seconds) / 60, 2), " min") # DWI Deep Learning Segmentation dwi_mask_sagittal = predict_mask(cases_file_s, trained_model_folder, view='sagittal') @@ -741,7 +729,7 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): end_masking_time = datetime.datetime.now() total_masking_time = end_masking_time - start_total_time - total_preprocessing_time - print ("Masking Time Taken : ", round(int(total_masking_time.seconds)/60, 2), " min") + print("Masking Time Taken : ", round(int(total_masking_time.seconds) / 60, 2), " min") transformed_file = registered_file omat_file = mat_file @@ -751,37 +739,35 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): omat_list = [line.rstrip('\n') for line in open(omat_file)] # Post Processing - print ("Splitting files....") + print("Splitting files....") cases_mask_sagittal = split(dwi_mask_sagittal, target_list, view='sagittal') cases_mask_coronal = split(dwi_mask_coronal, target_list, view='coronal') cases_mask_axial = split(dwi_mask_axial, target_list, view='axial') multi_mask = [] for i in range(0, len(cases_mask_sagittal)): - sagittal_SO = cases_mask_sagittal[i] coronal_SO = cases_mask_coronal[i] axial_SO = cases_mask_axial[i] input_file = target_list[i] - multi_view_mask = multi_view_fast(sagittal_SO, - coronal_SO, - axial_SO, + multi_view_mask = multi_view_fast(sagittal_SO, + coronal_SO, + axial_SO, input_file) - - brain_mask_multi = npy_to_nifti(list(transformed_cases[i].split()), - list(multi_view_mask.split()), + brain_mask_multi = npy_to_nifti(list(transformed_cases[i].split()), + list(multi_view_mask.split()), list(target_list[i].split()), - view='multi', - reference=list(target_list[i].split()), + view='multi', + reference=list(target_list[i].split()), omat=list(omat_list[i].split())) - - print ("Mask file : ", brain_mask_multi) + print("Mask file : ", brain_mask_multi) multi_mask.append(brain_mask_multi[0]) - quality_control(multi_mask, target_list, tmp_path, view='multi') + if args.snap: + quality_control(multi_mask, target_list, tmp_path, view='multi') if args.Sagittal: omat = omat_list @@ -789,14 +775,15 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): omat = None if args.Sagittal: - sagittal_mask = npy_to_nifti(transformed_cases, - cases_mask_sagittal, - target_list, - view='sagittal', - reference=target_list, - omat=omat) + sagittal_mask = npy_to_nifti(transformed_cases, + cases_mask_sagittal, + target_list, + view='sagittal', + reference=target_list, + omat=omat) list_masks(sagittal_mask, view='sagittal') - quality_control(sagittal_mask, target_list, tmp_path, view='sagittal') + if args.snap: + quality_control(sagittal_mask, target_list, tmp_path, view='sagittal') if args.Coronal: omat = omat_list @@ -804,14 +791,15 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): omat = None if args.Coronal: - coronal_mask = npy_to_nifti(transformed_cases, - cases_mask_coronal, - target_list, - view='coronal', - reference=target_list, - omat=omat) + coronal_mask = npy_to_nifti(transformed_cases, + cases_mask_coronal, + target_list, + view='coronal', + reference=target_list, + omat=omat) list_masks(coronal_mask, view='coronal') - quality_control(coronal_mask, target_list, tmp_path, view='coronal') + if args.snap: + quality_control(coronal_mask, target_list, tmp_path, view='coronal') if args.Axial: omat = omat_list @@ -819,14 +807,15 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): omat = None if args.Axial: - axial_mask = npy_to_nifti(transformed_cases, - cases_mask_axial, - target_list, - view='axial', - reference=target_list, - omat=omat) + axial_mask = npy_to_nifti(transformed_cases, + cases_mask_axial, + target_list, + view='axial', + reference=target_list, + omat=omat) list_masks(axial_mask, view='axial') - quality_control(axial_mask, target_list, tmp_path, view='axial') + if args.snap: + quality_control(axial_mask, target_list, tmp_path, view='axial') for i in range(0, len(cases_mask_sagittal)): clear(path.dirname(cases_mask_sagittal[i])) @@ -842,5 +831,4 @@ def quality_control(mask_list, target_list, tmp_path, view='default'): end_total_time = datetime.datetime.now() total_t = end_total_time - start_total_time - print ("Total Time Taken : ", round(int(total_t.seconds)/60, 2), " min") - + print("Total Time Taken : ", round(int(total_t.seconds) / 60, 2), " min") diff --git a/src/maskfilter.py b/src/maskfilter.py index 8ee656e..7e0bd00 100755 --- a/src/maskfilter.py +++ b/src/maskfilter.py @@ -79,7 +79,7 @@ def maskfilter(maskPath, scale, filtered_maskPath): mask= nib.load(maskPath) - filtered_mask= perform_morph(mask.get_data(), scale) + filtered_mask= perform_morph(mask.get_fdata(), scale) nib.Nifti1Image(filtered_mask, affine= mask.affine, header= mask.header).to_filename(filtered_maskPath) diff --git a/src/preprocess_b0.py b/src/preprocess_b0.py index 84f242f..ceddd76 100755 --- a/src/preprocess_b0.py +++ b/src/preprocess_b0.py @@ -12,7 +12,7 @@ def process_trainingdata(dwib0_arr): count = 0 for b0 in dwib0_arr: img = nib.load(b0) - imgF32 = img.get_data().astype(np.float32) + imgF32 = img.get_fdata().astype(np.float32) ''' Intensity based segmentation of MR images is hampered by radio frerquency field inhomogeneity causing intensity variation. The intensity range is typically scaled between the highest and lowest signal in the Image. Intensity values diff --git a/src/preprocess_mask.py b/src/preprocess_mask.py index 589be7f..e77e4f3 100755 --- a/src/preprocess_mask.py +++ b/src/preprocess_mask.py @@ -11,7 +11,7 @@ def process_trainingdata(mask_arr): count = 0 for b0_mask in mask_arr: img = nib.load(b0_mask) - imgU8_sagittal = img.get_data().astype(np.uint8) # sagittal view + imgU8_sagittal = img.get_fdata().astype(np.uint8) # sagittal view imgU8_sagittal[ imgU8_sagittal < 0 ] = 0 imgU8_sagittal[ imgU8_sagittal > 1 ] = 1 imgU8_coronal = np.swapaxes(imgU8_sagittal,0,1) # coronal view diff --git a/tests/README.md b/tests/README.md new file mode 100644 index 0000000..8697d6b --- /dev/null +++ b/tests/README.md @@ -0,0 +1,22 @@ +### Tests + +[pipeline_test.sh](./pipeline_test.sh) is a convenient testing infrastructure. +In addition to the prerequisites mentioned in [instructions](../docs/README.md#installation), +it is depdendent on the following executables: + +* dcm2niix +* fslroi +* bet + +We recommend only putting the above executables in `PATH` without sourcing extenal environments +such as FSL's. For example--`fslroi`, `bet` could be put in `PATH` as: + +> export PATH=/path/to/fsl-6.0.7/share/fsl/bin:$PATH + +In addition, you need to set: + +> export FSLOUTPUTTYPE=NIFTI_GZ + +Moreover, you can put an `exit` statement after any line between 70-81 to run one or more tests: +https://github.com/pnlbwh/CNN-Diffusion-MRIBrain-Segmentation/blob/a81ef7e939714f88b67c0a6e84a0ff6db7004622/tests/pipeline_test.sh#L70-L81 +