8-3-UNet提取海马体

字数统计: 4.3k字   |   阅读时长: 22分
× 文章目录
  1. 1. 1-横截面选取
  2. 2. 2-同组内混合
  3. 3. 3-横截面 Axial进行方位校正，以便海马体的提取
  4. 4. 4-UNet进行海马体分割
    1. 4.1. 搭建UNet
    2. 4.2. 模型对ADNI-Datasets 方向校正后的横截面数据作海马体预测提取
  5. 5. 5-海马体候选挑选
    1. 5.1. 挑选颜色均值>2的L_R图片
    2. 5.2. 融合图片,将左右海马体融合到原图上
  6. 6. 6-补充AD-MCI
    1. 6.1. FSL左右海马体分割+freesurfer合并
      1. 6.1.1. 背景
      2. 6.1.2. bash直接运行
      3. 6.1.3. python+nipype实现

Author: haoransun
Wechat: SHR—97

1-横截面选取

Axial：又名横断面，transverseCoronal：又名冠状面，Sagittal：矢状面，如同一个箭矢劈开成左右两半

image

这种图就是横断面Axial

image

nii_process-origin.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
import scipy, numpy, shutil, os, nibabel, glob
import sys, getopt
import imageio

# 切片 矢状面进行切片（start=80, end=110 total 30）
def niito2D_A_batch(nii_file, slice_dim, parent_folder):
    
    outputfile = os.path.join(nii_file.split("/")[-1][: 15], slice_dim)
    outputfile = os.path.join(parent_folder, outputfile)       #输出文件夹

    image_array = nibabel.load(nii_file).get_data() #数据读取

    # set destination folder
    if not os.path.exists(outputfile):
        os.makedirs(outputfile)   #不存在输出文件夹则新建

    slice_start = 80 #从第几个切片开始
    slice_end = 110

    # iterate through slices
    for current_slice in range(slice_start, slice_end):
        # alternate slices
        # [A,B,C]: A 矢状面 B 冠状面 C 横断面
                
        data = image_array[:, :, current_slice] 
        #切片命名 
        # /home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin/ADNI_002_S_0619_I120964.nii
        # /home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin/ADNI_002_S_061081.png
        image_name = nii_file[:-12] +"_" "{:0>3}".format(str(current_slice + 1)) + ".png"

        #保存
        imageio.imwrite(image_name, data)
        src = image_name
        shutil.move(src, outputfile)


# 切片 冠状面进行切片（start=80, end=120 total 40）
def niito2D_B_batch(nii_file, slice_dim, parent_folder):
    
    outputfile = os.path.join(nii_file.split("/")[-1][: 15], slice_dim)
    outputfile = os.path.join(parent_folder, outputfile)       #输出文件夹

    
    image_array = nibabel.load(nii_file).get_data() #数据读取

    # set destination folder
    if not os.path.exists(outputfile):
        os.makedirs(outputfile)   #不存在输出文件夹则新建

    slice_start = 80 #从第几个切片开始
    slice_end = 120

    # iterate through slices
    for current_slice in range(slice_start, slice_end):
        # alternate slices
        # [A,B,C]: A 矢状面 B 冠状面 C 横断面
                
        data = image_array[:, current_slice, :] 
        #切片命名 
        # /home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin/ADNI_002_S_0619_I120964.nii
        # /home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin/ADNI_002_S_061081.png
        image_name = nii_file[:-12] +"_" "{:0>3}".format(str(current_slice + 1)) + ".png"

        #保存
        imageio.imwrite(image_name, data)
        src = image_name
        shutil.move(src, outputfile)


# 切片 横断面进行切片（start=100, end=130 total 30）
def niito2D_C_batch(nii_file, slice_dim, parent_folder):
    
    outputfile = os.path.join(nii_file.split("/")[-1][: 15], slice_dim)
    outputfile = os.path.join(parent_folder, outputfile)       #输出文件夹
    
    image_array = nibabel.load(nii_file).get_data() #数据读取

    # set destination folder
    if not os.path.exists(outputfile):
        os.makedirs(outputfile)   #不存在输出文件夹则新建

    slice_start = 100 #从第几个切片开始
    slice_end = 130

    # iterate through slices
    for current_slice in range(slice_start, slice_end):
        # alternate slices
        # [A,B,C]: A 矢状面 B 冠状面 C 横断面
                
        data = image_array[current_slice, :, :] 
        #切片命名 
        # /home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin/ADNI_002_S_0619_I120964.nii
        # /home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin/ADNI_002_S_061081.png
        image_name = nii_file[:-12] +"_" "{:0>3}".format(str(current_slice + 1)) + ".png"

        #保存
        imageio.imwrite(image_name, data)
        src = image_name
        shutil.move(src, outputfile)


AD_Brain_Affine_list = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin/*.nii")

for i in AD_Brain_Affine_list:
    # niito2D_A_batch(i, "slice_A", "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices/")
    # niito2D_B_batch(i, "slice_B", "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices/")
    niito2D_C_batch(i, "slice_C", "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices/")


CN_Brain_Affine_list = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin/*.nii")

for i in CN_Brain_Affine_list:
    # niito2D_A_batch(i, "slice_A", "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices/")
    # niito2D_B_batch(i, "slice_B", "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices/")
    niito2D_C_batch(i, "slice_C", "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices/")


MCI_Brain_Affine_list = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin/*.nii")

for i in MCI_Brain_Affine_list:
    # niito2D_A_batch(i, "slice_A", "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices/")
    # niito2D_B_batch(i, "slice_B", "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices/")
    niito2D_C_batch(i, "slice_C", "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices/")

2-同组内混合

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
slices_AD_C = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices/*/slice_C/*.png")
target_AD_C = "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C/"

for i in slices_AD_C:
    # shutil.copy(i, target_AD_C)

print(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C/*.png")))


slices_CN_C = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices/*/slice_C/*.png")
target_CN_C = "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C/"

for i in slices_CN_C:
    # shutil.copy(i, target_CN_C)

print(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C/*.png")))



slices_MCI_C = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices/*/slice_C/*.png")
target_MCI_C = "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C/"

for i in slices_MCI_C:
    # shutil.copy(i, target_MCI_C)

print(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C/*.png")))

3-横截面 Axial进行方位校正，以便海马体的提取

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
import pandas as pd
import numpy as np
import os
import shutil
import cv2
import glob
import imageio


slices_AD_C = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C/*.png")
target_AD_C = "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD/"
for i in slices_AD_C:
    img = cv2.imread(i)
    img = cv2.rotate(img, cv2.ROTATE_180)
    img_name = i.split("/")[-1]
    imageio.imwrite(img_name, img)
    shutil.move(img_name, target_AD_C)

print(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD/*.png")))
    
slices_CN_C = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C/*.png")
target_CN_C = "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD/"
for i in slices_CN_C:
    img = cv2.imread(i)
    img = cv2.rotate(img, cv2.ROTATE_180)
    img_name = i.split("/")[-1]
    imageio.imwrite(img_name, img)
    shutil.move(img_name, target_CN_C)

print(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD/*.png")))


slices_MCI_C = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C/*.png")
target_MCI_C = "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD/"
for i in slices_MCI_C:
    img = cv2.imread(i)
    img = cv2.rotate(img, cv2.ROTATE_180)
    img_name = i.split("/")[-1]
    imageio.imwrite(img_name, img)
    shutil.move(img_name, target_MCI_C)

print(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD/*.png")))

4-UNet进行海马体分割

搭建UNet

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import Dataset
from torch.utils.data import DataLoader
from torch.utils.data.dataset import random_split
from torchvision import transforms
import cv2
import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
import os
import glob

class MRIDataset(Dataset):
    def __init__(self, url, transform=None):
        self.imgs = glob.glob(os.path.join(url, "img/*.jpg"))
        self.transform = transform
    
    def augment(self, img, flip):
        # flip = 1 水平翻转
        # flip = 2 垂直翻转
        # flip = 3 同时进行水平翻转和垂直翻转
        img_flipped = cv2.flip(img, flip)
        return img_flipped
    
    def __getitem__(self, index):
        
        # 通过索引获取原图和标签的URL
        img_url = self.imgs[index]
        label_left_url = img_url.replace("img", "label").replace("ACPC", "L")
        label_right_url = img_url.replace("img", "label").replace("ACPC", "R")
        
        # 通过cv2库读取图像
        image = cv2.imread(img_url)
        label_left = cv2.imread(label_left_url)
        label_right = cv2.imread(label_right_url)
        
                
        width = 224
        height = 224
        image = cv2.resize(image,(width,height))
        label_left = cv2.resize(label_left,(width,height))
        label_right = cv2.resize(label_right,(width,height))
        
        
        bImg, gImg, rImg = cv2.split(image)
        image = cv2.merge([bImg,gImg,rImg])
        
        blabel_left, glabel_left, rlabel_left = cv2.split(label_left)
        label_left = cv2.merge([blabel_left,glabel_left,rlabel_left])
        
        blabel_right, glabel_right, rlabel_right = cv2.split(label_right)
        label_right = cv2.merge([blabel_right,glabel_right,rlabel_right])
        
        # 合并左右海马体
        label = cv2.addWeighted(label_left,1,label_right,1,0.2) 
        
        # 定位到转灰度时图片要为三通道，否则报错，在转灰度前，加上如下代码将图转为三通道即可
        # 采用CV_BGR2GRAY,  转换公式Gray = 0.1140B + 0.5870G + 0.2989*R
        # image = cv2.merge([image,image,image])
        # label = cv2.merge([label,label,label])
        
        # 转换为灰度图：3通道转换为1通道
        image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
        label = cv2.cvtColor(label, cv2.COLOR_BGR2GRAY)
        
        # 修正标签，将海马体部分像素设置为1，非海马体部分设置为0
        if label.max()>1:
            label = label/255
        label[label>=0.5] = 1
        label[label<0.5] = 0
        
        # 进行随机的数据增强
        flip = random.choice([-1, 0, 1, 2])
        if flip != 2:
            image = self.augment(image, flip)
            label = self.augment(label, flip)
        
        # 转换数据
        if self.transform:
            image = self.transform(image)
            label = self.transform(label)

        return image, label
    
    def __len__(self):
        return len(self.imgs)


class ConvUnit(nn.Module):
    def __init__(self, in_channels, out_channels):
        super(ConvUnit, self).__init__()
        self.unit = nn.Sequential(
            # 保持图像大小
            nn.Conv2d(in_channels, out_channels, kernel_size=3, stride=1, padding=1),
            # 常用于卷积网络中卷积操作之后，进行数据的归一化处理(防止梯度消失或爆炸)
            # 使得数据在进行Relu之前不会因为数据过大而导致网络性能不稳定
            nn.BatchNorm2d(out_channels),
            # inplace=True, 地址传递，会改变输入数据的值,节省反复申请与释放内存的空间与时间
            nn.ReLU(inplace=True),
            
             # 保持图像大小
            nn.Conv2d(out_channels, out_channels, kernel_size=3, stride=1, padding=1),
            nn.BatchNorm2d(out_channels),
            nn.ReLU(inplace=True)
        )
        
    def forward(self, x):
        return self.unit(x)


class DownSampling(nn.Module):
    def __init__(self, in_channels, out_channels):
        super(DownSampling, self).__init__()
        self.layer = nn.Sequential(
            # 最大池化，将图像大小变为原来的 1/2
            nn.MaxPool2d(2),
            # 卷积单元进行卷积操作
            ConvUnit(in_channels, out_channels)
        )
    
    def forward(self, x):
        return self.layer(x)


class UpSampling(nn.Module):
    def __init__(self, in_channels, out_channels):
        super(UpSampling, self).__init__()
        
        # 上采样层，将通道数变为1/2是为了在concat后保持通道数不变
        self.layer = nn.ConvTranspose2d(in_channels, in_channels//2, kernel_size=2, stride=2)
       
        # 卷积单元
        self.conv = ConvUnit(in_channels, out_channels)
    
    def forward(self, x, r):
        # 对x进行上采样，同时通道数减半
        x = self.layer(x)
        
        # 将x与r在通道维度连接，恢复原本通道数
        x = torch.cat((x, r), dim=1)
        
        return self.conv(x)


class UNet(nn.Module):
    def __init__(self, in_channels, out_channels):
        super(UNet, self).__init__()
        # 输入层
        self.conv = ConvUnit(in_channels, 64)
        # 定义四个下采样层和四个上采样层
        self.D1 = DownSampling(64, 128)
        self.D2 = DownSampling(128, 256)
        self.D3 = DownSampling(256, 512)
        self.D4 = DownSampling(512, 1024)
        self.U1 = UpSampling(1024, 512)
        self.U2 = UpSampling(512, 256)
        self.U3 = UpSampling(256, 128)
        self.U4 = UpSampling(128, 64)
        # 输出层，输出图像像素保持在0～1以进行二分类
        self.out = nn.Sequential(
            nn.Conv2d(64, out_channels, kernel_size=1),
            nn.Sigmoid()
        )
        
    def forward(self, x):
        # U-NET的左半部分
        L1 = self.conv(x)
        L2 = self.D1(L1)
        L3 = self.D2(L2)
        L4 = self.D3(L3)
        Bottom = self.D4(L4)
        # U-NET的右半部分
        R1 = self.U1(Bottom, L4)
        R2 = self.U2(R1, L3)
        R3 = self.U3(R2, L2)
        R4 = self.U4(R3, L1)
        return self.out(R4)

模型对ADNI-Datasets 方向校正后的横截面数据作海马体预测提取

1
2
3
4
5
6
7
8
9
import glob

print("AD: " + str(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD/*.png"))))
print("CN: " + str(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD/*.png"))))
print("MCI: "+ str(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD/*.png"))))

print("AD: " + str(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs/*.png"))))
print("CN: " + str(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs/*.png"))))
print("MCI: "+ str(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs/*.png"))))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
class MRI_3_class_Dataset(Dataset):
    def __init__(self, url,  transform=None):
        self.imgs = glob.glob(url)
        self.transform = transform
    
    def augment(self, img, flip):
        # flip = 1 水平翻转
        # flip = 2 垂直翻转
        # flip = 3 同时进行水平翻转和垂直翻转
        img_flipped = cv2.flip(img, flip)
        return img_flipped
    
    def __getitem__(self, index):
        
        # 通过索引获取原图和标签的URL
        img_url = self.imgs[index]
       
        # 通过cv2库读取图像
        image = cv2.imread(img_url)
                
        width = 224
        height = 224
        image = cv2.resize(image,(width,height))
        
        bImg, gImg, rImg = cv2.split(image)
        image = cv2.merge([bImg,gImg,rImg])
        
        # 定位到转灰度时图片要为三通道，否则报错，在转灰度前，加上如下代码将图转为三通道即可
        # 采用CV_BGR2GRAY,  转换公式Gray = 0.1140B + 0.5870G + 0.2989*R
        # image = cv2.merge([image,image,image])
        # label = cv2.merge([label,label,label])
        
        # 转换为灰度图：3通道转换为1通道
        image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
        
        # 进行随机的数据增强
        # flip = random.choice([-1, 0, 1, 2])
        # if flip != 2:
        #     image = self.augment(image, flip)
        
        # 转换数据
        if self.transform:
            image = self.transform(image)
            
        return image
    
    def __len__(self):
        return len(self.imgs)



def load_mri_3_class_data(url, batch_size=4, shuffle=True, split=None):
    
    # 实例化自定义Dataset，加载数据集
    mri = MRI_3_class_Dataset(url, transform=transforms.Compose([
        # ToTensor()将shape为(H, W, C)的nump.ndarray或img转为shape为(C, H, W)的tensor，其将每一个数值归一化到[0,1]
        # 其归一化方法比较简单，直接除以255即可
        transforms.ToTensor()
    ]))
    
    # 将数据集封装到DataLoader
    data_loader = DataLoader(mri, batch_size=batch_size, shuffle=shuffle)
    return data_loader

预训练好的模型在 阿里云 Unet中

1
model = "/home/pugongying/data/jupyterLab/shr/AD/U-NET/model_version_2/model.pth"
1
2
3
4
5
6
7
8
9
10
AD_Data_Loader = load_mri_3_class_data(url="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD/*.png",   batch_size=4, shuffle=True, split=None)
predict(AD_Data_Loader, model=model, pred_dir= "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs")


CN_Data_Loader = load_mri_3_class_data(url="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD/*.png",   batch_size=4, shuffle=True, split=None)
predict(CN_Data_Loader, model=model, pred_dir= "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs")


MCI_Data_Loader = load_mri_3_class_data(url="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD/*.png",   batch_size=4, shuffle=True, split=None)
predict(MCI_Data_Loader, model=model, pred_dir= "/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs")

5-海马体候选挑选

1
2
进入命令行， cp -r origin-slices_C_CD_hs/*.png      origin-slices_C_CD_hs_
因为使用Unet分割一次太慢，所以需要保留分割结果

挑选颜色均值>2的L_R图片

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
import os
import glob
import shutil
import numpy as np
import pandas as pd
import random
import cv2
import matplotlib.pyplot as plt
from PIL import Image
import skimage




AD_Origin_LR_HS= glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs/*_L_R.png")
CN_Origin_LR_HS= glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs/*_L_R.png")
MCI_Origin_LR_HS= glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs/*_L_R.png")

print("AD_Origin_LR_HS: " + str(len(AD_Origin_LR_HS)))
print("CN_Origin_LR_HS: " + str(len(CN_Origin_LR_HS)))
print("MCI_Origin_LR_HS: " + str(len(MCI_Origin_LR_HS)))
print()


AD_Origin_LR_HS_= glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs_/*_L_R.png")
CN_Origin_LR_HS_= glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs_/*_L_R.png")
MCI_Origin_LR_HS_= glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs_/*_L_R.png")

print("AD_Origin_LR_HS_: " + str(len(AD_Origin_LR_HS_)))
print("CN_Origin_LR_HS_: " + str(len(CN_Origin_LR_HS_)))
print("MCI_Origin_LR_HS_: " + str(len(MCI_Origin_LR_HS_)))
print()


AD_Origin_HS= glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs/*.png")
CN_Origin_HS= glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs/*.png")
MCI_Origin_HS= glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs/*.png")


print("AD_Origin_HS: " + str(len(AD_Origin_HS)))
print("CN_Origin_HS: " + str(len(CN_Origin_HS)))
print("MCI_Origin_HS: " + str(len(MCI_Origin_HS)))




# 遍历 AD-CN-MCI数据集，白点出现超过 才可以留存
def select_candidate_first(url):
    labels = glob.glob(url)
    for candidate_label in labels:
        # 对图片进行判断，图片平均值>1(实际平均值要为0，此处提供的图片全黑为1)，就是全黑，否则就是有白点
        origin_image = candidate_label.replace("_L_R", "")
        if (np.mean(cv2.imread(candidate_label))) >2:
            pass
        else:
            os.remove(candidate_label)
            os.remove(origin_image)




select_candidate_first(url="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs_/*_L_R.png")
select_candidate_first(url="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs_/*_L_R.png")
select_candidate_first(url="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs_/*_L_R.png")

test_image = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs_/*_L_R.png")
test_image.sort()

plt.figure(figsize=(32,32))

k = 1
for i in test_image:
    
    image_mean_score = np.mean(cv2.imread(i))
    print(str(i[-7:-4]) +" : " +  str(image_mean_score))
    
    plt.subplot(10,6,k)
    plt.title(image_mean_score)
    plt.imshow(Image.open(i))
    plt.legend()
    
    k = k+1



test_image = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs_/*_L_R.png")
test_image.sort()

plt.figure(figsize=(32,32))

k = 1
for i in test_image:
    
    image_mean_score = np.mean(cv2.imread(i))
    print(str(i[-7:-4]) +" : " +  str(image_mean_score))
    
    plt.subplot(10,6,k)
    plt.title(image_mean_score)
    plt.imshow(Image.open(i))
    plt.legend()
    
    k = k+1



test_image = glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs_/*_L_R.png")
test_image.sort()

plt.figure(figsize=(32,32))

k = 1
for i in test_image:
    
    image_mean_score = np.mean(cv2.imread(i))
    print(str(i[-7:-4]) +" : " +  str(image_mean_score))
    
    plt.subplot(10,6,k)
    plt.title(image_mean_score)
    plt.imshow(Image.open(i))
    plt.legend()
    
    k = k+1

融合图片,将左右海马体融合到原图上

  • 融合后原图放在origin-slices_C_CD_hs_Mix中
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def concat_iamge_hippocampus(url, target):
    labels = glob.glob(url)
    for candidate_label in labels:
        candidate_hs_image = cv2.imread(candidate_label)
        origin_img = candidate_label.replace("_L_R", "")
        origin_img_img = origin_img.split("/")[-1]
        img = cv2.imread(origin_img)
        main_concat = cv2.addWeighted(img[:,:,:], 1, candidate_hs_image[:,:,:], 1, 0.2)
        cv2.imwrite(os.path.join(target, origin_img_img), main_concat)


concat_iamge_hippocampus(url="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs_/*_L_R.png", target="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs_Mix")
print("AD_Origin_Mix: " + str(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/AD/origin-slices_C_CD_hs_Mix/*.png"))))


concat_iamge_hippocampus(url="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs_/*_L_R.png", target="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs_Mix")
print("CN_Origin_Mix: " + str(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/CN/origin-slices_C_CD_hs_Mix/*.png"))))



concat_iamge_hippocampus(url="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs_/*_L_R.png", target="/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs_Mix")
print("MCI_Origin_Mix: " + str(len(glob.glob("/home/pugongying/data/jupyterLab/shr/AD/ADNI-Datasets/MCI/origin-slices_C_CD_hs_Mix/*.png"))))

6-补充AD-MCI

AD_Origin_Mix: 178

CN_Origin_Mix: 352

MCI_Origin_Mix: 215

image

所以需要再多下载AD和MCI 进行补充

补充后：

AD_Origin_Mix: 300

CN_Origin_Mix: 500

MCI_Origin_Mix: 380

FSL左右海马体分割+freesurfer合并

https://blog.csdn.net/qq_38851184/article/details/124452709

背景

需要分割出海马体。可借鉴方法：深度学习（UNet分割），形态学上的开闭，fsl中的分割。
使用fsl种的分割时需要注意，fsl分割分为左海马和右海马
方法：使用fsl中的first命令进行分割。
先在终端输入first查看需要输入的参数：

image

可以看到-i 要从哪个文件进行提取， -k输出分割文件名称，-m 分割的标准文件（*.bmv），-l

做flirt时的变换矩阵。

关于*.bmv文件：fsl自带的，在fsl安装目录data/first下面：

image

关于这两个文件夹有什么区别不清楚，但是海马在选中的文件夹中。

bash直接运行

flirt的变换矩阵即subjmat.mat,该文件在做配准的时候是需要用参数进行输出的，但是直接使用命令是默认不输出

命令：first -i subj.nii -l subjmat.mat -m L_Hipp_bin.bmv -k l_hippo.nii

同样的方法生成右海马：

first -i subj.nii -l subjmat.mat -m R_Hipp_bin.bmv -k r_hippo.nii

使用freesurfer中的mri_concat命令将左海马和右海马合并成整个海马体：

mri_concat –combine l_hippo.nii r_hippo.nii –o hippo.nii

得到hippo.nii就是海马文件：

image

python+nipype实现

分割成左右海马体+合并代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import os
import subprocess
from PyQt5 import QtWidgets
from PyQt5.QtWidgets import QMessageBox

def hippoSeg(in_file, out_file, out_path):
    '''
	segmentation hippocampus, and show the messagebox.
    :param in_file: subj-T1.nii
    :param out_file: hippo.nii
    :param out_path: /home/xxx/Desktop/test
    :return:
    '''

    l_hippo = os.path.join(out_path, 'L_Hipp_bin.bmv')
    r_hippo = os.path.join(out_path, 'R_Hipp_bin.bmv')

    mat_file = os.path.join(out_path, 'reg.mat') # 配准时输出的矩阵
    l_hippo_out = os.path.join(out_path, 'l_hippo.nii.gz')
    r_hippo_out = os.path.join(out_path, 'r_hippo.nii.gz')

    try:
   		''' 分割左海马体 '''
        l_hippo_cmd = 'first -i {0} -l {1} -m {2} -k {3}'.format(in_file, mat_file, l_hippo, l_hippo_out)
        # 因为是连接的远程环境，所以每次使用fsl中的命令需要重新inport一下，感觉有点繁琐。
        l_cmdline = 'export FSLDIR=/etc/fsl/5.0 && source $FSLDIR/fsl.sh && {0}'.format(l_hippo_cmd)
        task = subprocess.run(l_cmdline, shell=True, stdout=subprocess.PIPE)
    except BaseException as err:
        print('err\n', err)
        QMessageBox.information(QtWidgets.QWidget(), 'messagebox', 'left hippocampus segmentation error, see log file',
                                QMessageBox.Yes)
        return -1

    try:
    	''' 分割右海马体 '''
        r_hippo_cmd = 'first -i {0} -l {1} -m {2} -k {3}'.format(in_file, mat_file, r_hippo, r_hippo_out)
        r_cmdline = 'export FSLDIR=/etc/fsl/5.0 && source $FSLDIR/fsl.sh && {0}'.format(r_hippo_cmd)
        task = subprocess.run(r_cmdline, shell=True, stdout=subprocess.PIPE)
    except BaseException as err:
        print('err\n', err)
        QMessageBox.information(QtWidgets.QWidget(), 'messagebox', 'right hippocampus segmentation error, see log file',
                                QMessageBox.Yes)
        return -1

    try:
    	''' 左右海马体联合 '''
    	# 这条命令可以换成nipype包中的API接口，这里偷懒没有换。
        concat_cmd = 'mri_concat --combine {0} {1} --o {2}'.format(os.path.join(out_path, l_hippo_out), 
                                                                   os.path.join(out_path, r_hippo_out),
                                                                   os.path.join(out_path, out_file))
         ''' 同样需要export一下freesurfer '''
        merge_cmdline = 'export FREESURFER_HOME=/usr/local/freesurfer && source $FREESURFER_HOME/SetUpFreeSurfer.sh &&  {0}'.format(concat_cmd)
        print('start merge')
        task = subprocess.run(merge_cmdline, shell=True, stdout=subprocess.PIPE)
        print('finished!')
        QMessageBox.information(QtWidgets.QWidget(), 'messagebox', 'finished hippocampus ',
                                QMessageBox.Yes)
        return 1
    except BaseException as err:
        print('err\n', err)
        QMessageBox.warning(QtWidgets.QWidget(), 'messagebox', 'hippocampus merge error, see log file',
                                QMessageBox.Yes)
        return -1
1
first -i ADNI_002_S_0619_affine.nii.gz -l subject.mat -m /usr/local/fsl/data/first/models_317_bin/L_Hipp_bin.bmv -k l_hippo.nii

扫一扫，分享到微信

微信分享二维码
本站总访问量次 | 本站访客数人

tag:

    缺失模块。
    1、请确保node版本大于6.2
    2、在博客根目录（注意不是yilia根目录）执行以下命令：
    npm i hexo-generator-json-content --save

    3、在根目录_config.yml里添加配置： 

      jsonContent:
    meta: false
    pages: false
    posts:
      title: true
      date: true
      path: true
      text: false
      raw: false
      content: false
      slug: false
      updated: false
      comments: false
      link: false
      permalink: false
      excerpt: false
      categories: false
      tags: true

一个脱离了高级趣味的人!