使用MONAI和Nilearn加载和处理fMRI数据：一个逐步教程-Python教程-PHP中文网

使用monai和nilearn加载和处理fmri数据：一个逐步教程

本文档旨在指导初学者如何使用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]这行代码用于从文件名中提取受试者名称，因此请确保文件名格式正确。

协和·太初

国内首个针对罕见病领域的AI大模型

查看详情

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个核心并行处理
    # ... (其他代码)

登录后复制

注意：