本文档旨在指导初学者如何使用MONAI框架结合Nilearn库加载、预处理和存储fMRI数据。通过详细的代码示例和解释，您将学习如何将NIfTI格式的fMRI数据集成到现有的MONAI处理流程中，并进行必要的预处理步骤，例如裁剪、标准化和格式转换，最终将处理后的数据保存为PyTorch张量。
1. 环境准备
在开始之前，请确保您已安装以下必要的Python库：
- MONAI: 用于医学图像分析的PyTorch基础框架。
- Nilearn: 用于神经影像分析的Python模块，提供加载和处理NIfTI图像的功能。
- PyTorch: 用于张量计算的深度学习框架。
- Joblib: 用于简化并行计算的库。
您可以使用pip安装这些库：
pip install monai nilearn torch joblib
2. 数据加载
原始代码使用MONAI的LoadImage转换来加载图像。为了更好地兼容NIfTI文件并利用Nilearn的优势，我们将使用Nilearn的load_img函数来加载数据。
from nilearn.image import load_img
import numpy as np
def load_fmri_data(file_path):
"""
使用Nilearn加载fMRI数据并返回NumPy数组。
"""
try:
nifti_image = load_img(file_path)
data = nifti_image.get_fdata()
return data
except Exception as e:
print(f"加载文件时出错: {file_path}, 错误信息: {e}")
return None此函数使用nilearn.image.load_img加载NIfTI文件，然后使用get_fdata()方法将其转换为NumPy数组。如果加载过程中发生错误，则会打印错误信息并返回None。
3. 集成到MONAI流程
接下来，我们将修改原始代码中的read_data函数，使用新的load_fmri_data函数加载数据。
from monai.transforms import LoadImage
import torch
import os
import time
from multiprocessing import Process, Queue
from nilearn.image import load_img
import numpy as np
def read_data(filename,load_root,save_root,subj_name,count,queue=None,scaling_method=None, fill_zeroback=False):
print("processing: " + filename, flush=True)
path = os.path.join(load_root, filename)
try:
# load each nifti file
# data, meta = LoadImage()(path) # Original code
data = load_fmri_data(path) # Modified code
if data is None:
return None
except Exception as e:
print(f"加载文件时出错: {path}, 错误信息: {e}")
return None
#change this line according to your file names
save_dir = os.path.join(save_root,subj_name)
isExist = os.path.exists(save_dir)
if not isExist:
os.makedirs(save_dir)
# change this line according to your dataset
data = data[:, 14:-7, :, :]
# width, height, depth, time
# Inspect the fMRI file first using your visualization tool.
# Limit the ranges of width, height, and depth to be under 96. Crop the background, not the brain regions.
# Each dimension of fMRI registered to MNI space (2mm) is expected to be around 100.
# You can do this when you load each volume at the Dataset class, including padding backgrounds to fill dimensions under 96.
background = data==0
if scaling_method == 'z-norm':
global_mean = data[~background].mean()
global_std = data[~background].std()
data_temp = (data - global_mean) / global_std
elif scaling_method == 'minmax':
data_temp = (data - data[~background].min()) / (data[~background].max() - data[~background].min())
data_global = torch.empty(data.shape)
data_global[background] = data_temp[~background].min() if not fill_zeroback else 0
# data_temp[~background].min() is expected to be 0 for scaling_method == 'minmax', and minimum z-value for scaling_method == 'z-norm'
data_global[~background] = data_temp[~background]
# save volumes one-by-one in fp16 format.
data_global = torch.tensor(data_global).type(torch.float16) #Convert numpy array to tensor.
data_global_split = torch.split(data_global, 1, 3)
for i, TR in enumerate(data_global_split):
torch.save(TR.clone(), os.path.join(save_dir,"frame_"+str(i)+".pt"))关键修改：
- 注释掉 data, meta = LoadImage()(path) 这一行，并使用 data = load_fmri_data(path) 代替。
- 在将data_global保存为tensor时，需要将numpy array 转换为tensor。data_global = torch.tensor(data_global).type(torch.float16)
4. 目录结构和文件名
确保您的fMRI数据存储在正确的目录结构中，并且文件名符合代码的预期。根据原始代码，数据应该位于load_root指定的目录下，每个受试者的文件名为Sub1.nii等。subj_name = filename[:-7]这行代码用于从文件名中提取受试者名称，因此请确保文件名格式正确。
5. 完整代码示例
以下是集成了Nilearn的完整代码示例：
from monai.transforms import LoadImage
import torch
import os
import time
from multiprocessing import Process, Queue
from nilearn.image import load_img
import numpy as np
def load_fmri_data(file_path):
"""
使用Nilearn加载fMRI数据并返回NumPy数组。
"""
try:
nifti_image = load_img(file_path)
data = nifti_image.get_fdata()
return data
except Exception as e:
print(f"加载文件时出错: {file_path}, 错误信息: {e}")
return None
def read_data(filename,load_root,save_root,subj_name,count,queue=None,scaling_method=None, fill_zeroback=False):
print("processing: " + filename, flush=True)
path = os.path.join(load_root, filename)
try:
# load each nifti file
# data, meta = LoadImage()(path) # Original code
data = load_fmri_data(path) # Modified code
if data is None:
return None
except Exception as e:
print(f"加载文件时出错: {path}, 错误信息: {e}")
return None
#change this line according to your file names
save_dir = os.path.join(save_root,subj_name)
isExist = os.path.exists(save_dir)
if not isExist:
os.makedirs(save_dir)
# change this line according to your dataset
data = data[:, 14:-7, :, :]
# width, height, depth, time
# Inspect the fMRI file first using your visualization tool.
# Limit the ranges of width, height, and depth to be under 96. Crop the background, not the brain regions.
# Each dimension of fMRI registered to MNI space (2mm) is expected to be around 100.
# You can do this when you load each volume at the Dataset class, including padding backgrounds to fill dimensions under 96.
background = data==0
if scaling_method == 'z-norm':
global_mean = data[~background].mean()
global_std = data[~background].std()
data_temp = (data - global_mean) / global_std
elif scaling_method == 'minmax':
data_temp = (data - data[~background].min()) / (data[~background].max() - data[~background].min())
data_global = np.empty(data.shape)
data_global[background] = data_temp[~background].min() if not fill_zeroback else 0
# data_temp[~background].min() is expected to be 0 for scaling_method == 'minmax', and minimum z-value for scaling_method == 'z-norm'
data_global[~background] = data_temp[~background]
# save volumes one-by-one in fp16 format.
data_global = torch.tensor(data_global).type(torch.float16) #Convert numpy array to tensor.
data_global_split = torch.split(data_global, 1, 3)
for i, TR in enumerate(data_global_split):
torch.save(TR.clone(), os.path.join(save_dir,"frame_"+str(i)+".pt"))
def main():
# change two lines below according to your dataset
dataset_name = 'ABCD'
load_root = '/storage/4.cleaned_image' # This folder should have fMRI files in nifti format with subject names. Ex) sub-01.nii.gz
save_root = f'/storage/7.{dataset_name}_MNI_to_TRs_minmax'
scaling_method = 'z-norm' # choose either 'z-norm'(default) or 'minmax'.
# make result folders
filenames = os.listdir(load_root)
os.makedirs(os.path.join(save_root,'img'), exist_ok = True)
os.makedirs(os.path.join(save_root,'metadata'), exist_ok = True) # locate your metadata file at this folder
save_root = os.path.join(save_root,'img')
finished_samples = os.listdir(save_root)
queue = Queue()
count = 0
for filename in sorted(filenames):
subj_name = filename[:-4] # change the line to remove .nii
# extract subject name from nifti file. [:-7] rules out '.nii.gz'
# we recommend you use subj_name that aligns with the subject key in a metadata file.
expected_seq_length = 1000 # Specify the expected sequence length of fMRI for the case your preprocessing stopped unexpectedly and you try to resume the preprocessing.
# change the line below according to your folder structure
if (subj_name not in finished_samples) or (len(os.listdir(os.path.join(save_root,subj_name))) < expected_seq_length): # preprocess if the subject folder does not exist, or the number of pth files is lower than expected sequence length.
try:
count+=1
p = Process(target=read_data, args=(filename,load_root,save_root,subj_name,count,queue,scaling_method))
p.start()
if count % 32 == 0: # requires more than 32 cpu cores for parallel processing
p.join()
except Exception:
print('encountered problem with'+filename)
print(Exception)
if __name__=='__main__':
start_time = time.time()
main()
end_time = time.time()
print('\nTotal', round((end_time - start_time) / 60), 'minutes elapsed.')注意：
- 根据您的实际数据路径修改load_root和save_root。
- 根据您的文件名格式修改subj_name = filename[:-7]这一行。 例如，如果您的文件名为 S1.nii，则应该修改为 subj_name = filename[:-4]。
- 根据您的数据集调整data = data[:, 14:-7, :, :]这一行，以裁剪到感兴趣的脑区。
- 根据需要调整缩放方法（scaling_method）和缺失值填充方法（fill_zeroback）。
6. 并行处理
原始代码使用了multiprocessing模块进行并行处理。如果您需要进一步优化性能，可以考虑使用joblib库，它提供了更简洁的API和更好的并行化控制。
例如，您可以将read_data函数包装在一个Parallel循环中：
from joblib import Parallel, delayed
def process_file(filename, load_root, save_root, scaling_method):
subj_name = filename[:-4] # 调整后缀移除方式
read_data(filename, load_root, save_root, subj_name, 0, None, scaling_method)
def main():
# ... (其他代码)
filenames = os.listdir(load_root)
Parallel(n_jobs=4)(delayed(process_file)(filename, load_root, save_root, scaling_method) for filename in filenames) # 使用4个核心并行处理
# ... (其他代码)注意：
- 根据您的CPU核心数调整n_jobs参数。
- 并行处理可能会增加内存消耗，请确保您的系统有足够的内存。
7. 总结
通过将Nilearn集成到MONAI流程中，您可以更轻松地加载和处理fMRI数据。本文档提供了一个基本的示例，您可以根据您的具体需求进行修改和扩展。记住要仔细检查您的数据路径、文件名格式和预处理参数，以确保代码能够正确运行。
