-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathutils.py
185 lines (137 loc) · 6.36 KB
/
utils.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
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
181
182
183
184
185
# aini_utils.py
#
#
# (c) Ziyu Li, Qiyuan Tian, Artificial Intelligence in Neuroimaging Software, 2022
import os
import copy
import numpy as np
import pandas as pd
import nibabel as nb
import tensorflow.keras.backend as K
def block_ind(mask, sz_block=64, sz_pad=0):
# find indices of smallest block that covers whole brain
tmp = np.nonzero(mask)
xind = tmp[0]
yind = tmp[1]
zind = tmp[2]
xmin = np.min(xind);
xmax = np.max(xind)
ymin = np.min(yind);
ymax = np.max(yind)
zmin = np.min(zind);
zmax = np.max(zind)
ind_brain = [xmin, xmax, ymin, ymax, zmin, zmax]
# calculate number of blocks along each dimension
xlen = xmax - xmin + 1
ylen = ymax - ymin + 1
zlen = zmax - zmin + 1
nx = int(np.ceil(xlen / sz_block)) + sz_pad
ny = int(np.ceil(ylen / sz_block)) + sz_pad
nz = int(np.ceil(zlen / sz_block)) + sz_pad
# determine starting and ending indices of each block
xstart = xmin
ystart = ymin
zstart = zmin
xend = xmax - sz_block + 1
yend = ymax - sz_block + 1
zend = zmax - sz_block + 1
xind_block = np.round(np.linspace(xstart, xend, nx))
yind_block = np.round(np.linspace(ystart, yend, ny))
zind_block = np.round(np.linspace(zstart, zend, nz))
ind_block = np.zeros([xind_block.shape[0] * yind_block.shape[0] * zind_block.shape[0], 6])
count = 0
for ii in np.arange(0, xind_block.shape[0]):
for jj in np.arange(0, yind_block.shape[0]):
for kk in np.arange(0, zind_block.shape[0]):
ind_block[count, :] = np.array(
[xind_block[ii], xind_block[ii] + sz_block - 1, yind_block[jj], yind_block[jj] + sz_block - 1,
zind_block[kk], zind_block[kk] + sz_block - 1])
count = count + 1
ind_block = ind_block.astype(int)
return ind_block, ind_brain
def denormalize_image(imgall, imgnormall, mask): # transform the predicted image intensities to the original range
imgresall_denorm = np.zeros(imgall.shape)
for jj in np.arange(imgall.shape[-1]):
img = imgall[:, :, :, jj: jj + 1]
imgres = imgnormall[:, :, :, jj: jj + 1]
img_mean = np.mean(img[mask > 0.5])
img_std = np.std(img[mask > 0.5])
imgres_norm = (imgres * img_std + img_mean) * mask
imgresall_denorm[:, :, :, jj: jj + 1] = imgres_norm
return imgresall_denorm
def normalize_image(imgall, imgresall, mask, norm_ch='all'): # standardize the image intensities
imgall_norm = copy.deepcopy(imgall)
imgresall_norm = copy.deepcopy(imgresall)
if norm_ch == 'all':
norm_ch = np.arange(imgall.shape[-1])
for jj in norm_ch:
img = imgall[:, :, :, jj: jj + 1]
imgres = imgresall[:, :, :, jj: jj + 1]
img_mean = np.mean(img[mask > 0.5])
img_std = np.std(img[mask > 0.5])
img_norm = (img - img_mean) / img_std * mask
imgres_norm = (imgres - img_mean) / img_std * mask
imgall_norm[:, :, :, jj: jj + 1] = img_norm
imgresall_norm[:, :, :, jj: jj + 1] = imgres_norm
return imgall_norm, imgresall_norm
def extract_block(data, inds): # extract blocks from whole brain volume data
xsz_block = inds[0, 1] - inds[0, 0] + 1
ysz_block = inds[0, 3] - inds[0, 2] + 1
zsz_block = inds[0, 5] - inds[0, 4] + 1
ch_block = data.shape[-1]
blocks = np.zeros((inds.shape[0], xsz_block, ysz_block, zsz_block, ch_block))
# block_num = 1
for ii in np.arange(inds.shape[0]):
inds_this = inds[ii, :]
blocks[ii, :, :, :, :] = data[inds_this[0]:inds_this[1] + 1, inds_this[2]:inds_this[3] + 1,
inds_this[4]:inds_this[5] + 1, :]
return blocks
def mean_squared_error_weighted(y_true, y_pred):
loss_weights = y_true[:, :, :, :, -1:]
y_true_weighted = y_true[:, :, :, :, :-1] * loss_weights
y_pred_weighted = y_pred[:, :, :, :, :-1] * loss_weights
return K.mean(K.square(y_pred_weighted - y_true_weighted), axis=-1)
def mean_absolute_error_weighted(y_true, y_pred):
loss_weights = y_true[:, :, :, :, -1:]
y_true_weighted = y_true[:, :, :, :, :-1] * loss_weights
y_pred_weighted = y_pred[:, :, :, :, :-1] * loss_weights
return K.mean(K.abs(y_pred_weighted - y_true_weighted), axis=-1)
def block2brain(blocks, inds0, mask, sz_crop=0): # turn blocks into brain volume
vol_brain = np.zeros([mask.shape[0], mask.shape[1], mask.shape[2], blocks.shape[-1]])
vol_count = np.zeros([mask.shape[0], mask.shape[1], mask.shape[2], blocks.shape[-1]])
inds = copy.deepcopy(inds0)
indmax = np.max(inds, axis=0)
indmin = np.min(inds, axis=0)
for tt in np.arange(inds.shape[0]):
inds_this = inds[tt, :]
inds0_this = np.array([0, blocks.shape[1], 0, blocks.shape[2], 0, blocks.shape[3]])
if sz_crop > 0:
for mm in np.arange(0, 6, 2):
if inds_this[mm] == indmin[mm]:
continue
else:
inds_this[mm] = inds_this[mm] + sz_crop
inds0_this[mm] = inds0_this[mm] + sz_crop
for mm in np.arange(1, 6, 2):
if inds_this[mm] == indmax[mm]:
continue
else:
inds_this[mm] = inds_this[mm] - sz_crop
inds0_this[mm] = inds0_this[mm] - sz_crop
print(vol_brain.shape, blocks.shape)
print(inds_this, inds0_this)
vol_brain[inds_this[0]:inds_this[1] + 1, inds_this[2]:inds_this[3] + 1, inds_this[4]:inds_this[5] + 1, :] = \
vol_brain[inds_this[0]:inds_this[1] + 1, inds_this[2]:inds_this[3] + 1, inds_this[4]:inds_this[5] + 1, :] + \
blocks[tt, inds0_this[0]:inds0_this[1], inds0_this[2]:inds0_this[3], inds0_this[4]:inds0_this[5], :]
vol_count[inds_this[0]:inds_this[1] + 1, inds_this[2]:inds_this[3] + 1, inds_this[4]:inds_this[5] + 1, :] = \
vol_count[inds_this[0]:inds_this[1] + 1, inds_this[2]:inds_this[3] + 1, inds_this[4]:inds_this[5] + 1,
:] + 1.
vol_count[vol_count < 0.5] = 1.
vol_brain = vol_brain / vol_count
vol_brain = vol_brain * mask
vol_count = vol_count * mask
return vol_brain, vol_count
def save_nii(fpNii, data, fpRef):
new_header = header = nb.load(fpRef).header.copy()
new_img = nb.nifti1.Nifti1Image(data, None, header=new_header)
nb.save(new_img, fpNii)