/
mimosa_data.R
261 lines (233 loc) · 11.6 KB
/
mimosa_data.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
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
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
#' @title MIMoSA Training Data Frame
#'
#' @description This function creates the training vectors from a single MRI study that has FLAIR, T1, T2, and PD volumes as well as binary masks of lesions. The function can create a tissue mask for the data (or the user can supply a brain mask), the candidate voxels for lesion segmentation, smoothed volumes, and coupling maps. The user may supply already normalized data if they wish to use an alternative normalization method.
#' @param brain_mask brain or tissue mask of class nifti
#' @param FLAIR volume of class nifti
#' @param T1 volume of class nifti
#' @param T2 volume of class nifti. If not available use NULL.
#' @param PD volume of class nifti. If not available use NULL.
#' @param tissue is a logical value that determines whether the brain mask is a full brain mask or tissue mask (excludes CSF), should be FALSE unless you provide the tissue mask as the brain_mask object
#' @param gold_standard gold standard lesion segmentation mask of class nifti
#' @param normalize is 'no' by default and will not perform any normalization on data. To normalize data specify 'Z' for z-score normalization or 'WS' for WhiteStripe normalization
#' @param cand_mask is NULL to use candidate mask procedure proposed with method or a nifti object to be used as the candidate mask
#' @param slices vector of desired slices to train on, if NULL then train over the entire brain mask
#' @param orientation string value telling which orientation the training slices are specified in, can take the values of "axial", "sagittal", or "coronal"
#' @param cores 1 numeric indicating the number of cores to be used (no more than 4 is useful for this software implementation)
#' @param verbose logical indicating printing diagnostic output
#' @export
#' @importFrom fslr fslerode fsl_smooth
#' @importFrom neurobase check_nifti datatyper niftiarr zscore_img
#' @importFrom parallel mclapply
#' @importFrom oasis voxel_selection correct_image_dim
#' @importFrom stats cov.wt qnorm
#' @importFrom utils combn
#' @importFrom WhiteStripe whitestripe whitestripe_norm
#' @return List of objects
# @examples \dontrun{
#
#}
mimosa_data <- function (brain_mask, FLAIR, T1, T2 = NULL, PD = NULL, tissue = FALSE,
gold_standard = NULL, normalize = 'no', cand_mask = NULL, slices = NULL,
orientation = c("axial", "coronal", "sagittal"), cores = 1, verbose = TRUE) {
if (verbose) {
message("# Checking File inputs")
}
check_nifti2 <- function(x) {
if (is.null(x)) {
return(NULL)
}
else {
return(neurobase::check_nifti(x))
}
}
# Verify inputs are NIFTI objects
FLAIR = check_nifti2(FLAIR)
T1 = check_nifti2(T1)
T2 = check_nifti2(T2)
PD = check_nifti2(PD)
gold_standard = check_nifti2(gold_standard)
correct_image_dim2 <- function(x) {
if (is.null(x)) {
return(NULL)
}
else {
return(oasis::correct_image_dim(x))
}
}
# Correct image dimension in case originals do not match
FLAIR = correct_image_dim2(FLAIR)
T1 = correct_image_dim2(T1)
T2 = correct_image_dim2(T2)
PD = correct_image_dim2(PD)
gold_standard = correct_image_dim2(gold_standard)
# Assignment or creation of tissue mask and verification mask is 0/1
if(tissue == TRUE){
brain_mask = neurobase::check_nifti(brain_mask)
brain_mask = brain_mask > 0
brain_mask = neurobase::datatyper(brain_mask, trybyte = TRUE)
brain_mask = oasis::correct_image_dim(brain_mask)
tissue_mask = brain_mask
} else if (tissue == FALSE){
if(verbose){
message("# Getting Tissue Mask")
}
ero_brain_mask = fslr::fslerode(brain_mask, kopts = "-kernel box 5x5x5",
retimg = TRUE)
tissue_mask = oasis::voxel_selection(flair = FLAIR, brain_mask = ero_brain_mask, cutoff = 0.15)
}
# Make a list of the images
mimosa_data = list(FLAIR = FLAIR, T1 = T1, T2 = T2,
PD = PD)
rm(list = c("FLAIR","T1", "T2", "PD", "brain_mask"))
# Remove null values from list in case T2 or PD is specified as NULL
nulls = sapply(mimosa_data, is.null)
mimosa_data = mimosa_data[!nulls]
# Normalize by z-score approach if normalize = "Z"
if (normalize == 'Z') {
if (verbose) {
message("# Normalizing Images using Z-score")
}
mimosa_data = lapply(mimosa_data, zscore_img, mask = tissue_mask,
margin = NULL)
normalized = mimosa_data
nulls = sapply(normalized, is.null)
normalized = normalized[!nulls]
}
if (normalize == 'WS'){
if (verbose) {
message("# Normalizing Images using WhiteStripe")
}
# Run WhiteStripe on T1 and FLAIR
WS_T1 = WhiteStripe::whitestripe(mimosa_data$T1, type="T1", stripped = TRUE, verbose = verbose)
WS_T2 = WhiteStripe::whitestripe(mimosa_data$FLAIR, type="T2", stripped = TRUE, verose = verbose)
mimosa_data$T1 = WhiteStripe::whitestripe_norm(mimosa_data$T1, WS_T1$whitestripe.ind)
mimosa_data$FLAIR = WhiteStripe::whitestripe_norm(mimosa_data$FLAIR, WS_T2$whitestripe.ind)
if(is.null(T2)==FALSE & is.null(PD)==TRUE){
# Images include T1, FLAIR, T2
mimosa_data$T2 = WhiteStripe::whitestripe_norm(mimosa_data$T2, WS_T2$whitestripe.ind)
}
if(is.null(T2)==TRUE & is.null(PD)==FALSE){
# Images include T1, FLAIR, PD
mimosa_data$PD = WhiteStripe::whitestripe_norm(mimosa_data$PD, WS_T2$whitestripe.ind)
}
if(is.null(T2)==FALSE & is.null(PD)==FALSE){
# Images include T1, FLAIR, T2, PD
mimosa_data$T2 = WhiteStripe::whitestripe_norm(mimosa_data$T2, WS_T2$whitestripe.ind)
mimosa_data$PD = WhiteStripe::whitestripe_norm(mimosa_data$PD, WS_T2$whitestripe.ind)
}
}
# Obtain candidate voxels and create candidate mask
if (verbose) {
message("# Voxel Selection Procedure")
}
if(is.null(cand_mask)){
top_voxels = oasis::voxel_selection(flair = mimosa_data$FLAIR,
brain_mask = tissue_mask, cutoff = 0.85)
} else if(!is.null(cand_mask)){
cand_mask = neurobase::check_nifti(cand_mask)
top_voxels = cand_mask
}
if (verbose) {
message("# Smoothing Images: Sigma = 10")
}
# Obtain smoothed at 10 images
smoothed_mask = fslr::fsl_smooth(file = tissue_mask, sigma = 10,
smooth_mask = FALSE)
smooth_10 = mclapply(mimosa_data, fslr::fslsmooth, sigma = 10, mask = tissue_mask,
retimg = TRUE, smooth_mask = TRUE, smoothed_mask = smoothed_mask, mc.cores = cores)
if (verbose) {
message("# Smoothing Images: Sigma = 20")
}
# Obtain smoothed at 20 images
smoothed_mask = fslr::fsl_smooth(file = tissue_mask, sigma = 20,
smooth_mask = FALSE)
smooth_20 = mclapply(mimosa_data, fslr::fslsmooth, sigma = 20, mask = tissue_mask,
retimg = TRUE, smooth_mask = TRUE, smoothed_mask = smoothed_mask, mc.cores = cores)
# Add smoothed images to the mimosa_data list with proper names
img_names = names(mimosa_data)
names(smooth_10) = paste0(img_names, "_10")
mimosa_data = c(mimosa_data, smooth_10)
names(smooth_20) = paste0(img_names, "_20")
mimosa_data = c(mimosa_data, smooth_20)
# Remove the smoothed objects to conserve memory
smoothed = list(smooth_10 = smooth_10, smooth_20 = smooth_20)
rm(list = c("smooth_10", "smooth_20"))
# Initialize empty lists to store all the combinations of coupling intercepts and slopes
coupling_intercepts = list()
coupling_slopes = list()
combos = combn(img_names, 2)
list_names = as.data.frame(matrix(nrow = 1, ncol = dim(combos)[2]*2))
# Run coupling on all possible combinations of images
if (verbose) {
message("# Running Coupling")
}
for(i in 1:dim(combos)[2]){
temp_files = list(eval(parse(text = paste0('mimosa_data$', combos[1,i]))),
eval(parse(text = paste0('mimosa_data$', combos[2,i]))))
temp_return = imco(files=temp_files, brainMask=tissue_mask, subMask=top_voxels,
type="regression",
ref=1, fwhm=3, reverse=TRUE, verbose=verbose,
retimg=TRUE, outDir=NULL)
# Regressed Y on X and X on Y so create variable names to match: YonX_int, YonX_slope, XonY_int, XonY_slope
list_names[i] = paste0(combos[1,i], 'on' , combos[2,i])
list_names[(i+dim(combos)[2])] = paste0(combos[2,i], 'on' , combos[1,i])
coupling_intercepts[[i]] = temp_return$beta0
coupling_slopes[[i]] = temp_return$beta1
coupling_intercepts[[(i+dim(combos)[2])]] = temp_return$beta0_reverse
coupling_slopes[[(i+dim(combos)[2])]] = temp_return$beta1_reverse
if (verbose) {
message(paste0('# Ran Coupling for ', combos[1,i], ' and ' , combos[2,i], ' Combinations Successfully'))
}
}
# Rename coupling lists before appending with mimosa_data
names(coupling_intercepts) = paste0(list_names, '_intercepts')
names(coupling_slopes) = paste0(list_names, '_slopes')
# Append data
mimosa_data = append(mimosa_data, coupling_intercepts)
mimosa_data = append(mimosa_data, coupling_slopes)
mimosa_data$gold_standard = gold_standard
mimosa_data$top_voxels = top_voxels
# Transform full list of images into a dataframe
mimosa_dataframe = lapply(mimosa_data, function(X) masked = c(X[top_voxels == 1]))
mimosa_dataframe = do.call(cbind, mimosa_dataframe)
rownames(mimosa_dataframe) = NULL
# Bind indices of the top voxels to mimoda_data for reference
inds = neurobase::niftiarr(top_voxels, 1)*top_voxels
inds = which(inds > 0, arr.ind = TRUE)
colnames(inds) = c("axial", "coronal", "sagittal")
mimosa_dataframe = as.data.frame(cbind(inds, mimosa_dataframe))
mimosa_dataframe$top_voxels = NULL
# Allow for slices of images
if (!is.null(slices)) {
orientation = match.arg(orientation)
mimosa_dataframe = mimosa_dataframe[mimosa_dataframe[, orientation] %in% slices, ]
}
# Determine return object based on inputs
if(normalize == 'no' & tissue == TRUE){
# If normalize = FALSE we do not normalize, if tissue = true we treat the brain mask as the tissue_mask
## In this case do not return normalized images or the tissue mask as user already has them
return(list(mimosa_dataframe = mimosa_dataframe, top_voxels = top_voxels, smoothed = smoothed,
coupling_intercepts = coupling_intercepts, coupling_slopes = coupling_slopes))
}
if(normalize != 'no' & tissue == TRUE){
# If normalize a value then we normalize the images provided, if tissue is true we treat the brain mask as
## The tissue mask in this case return the normalized images but they have the tissue mask so do not return
return(list(mimosa_dataframe = mimosa_dataframe, top_voxels = top_voxels, smoothed = smoothed,
coupling_intercepts = coupling_intercepts, coupling_slopes = coupling_slopes,
normalized = normalized))
}
if(normalize == 'no' & tissue == FALSE){
# If normalize is FALSE then we normalize images if tissue is false we find the tissue mask
## Return only tissue
return(list(mimosa_dataframe = mimosa_dataframe, top_voxels = top_voxels, smoothed = smoothed,
coupling_intercepts = coupling_intercepts, coupling_slopes = coupling_slopes,
tissue_mask = tissue_mask))
}
if(normalize != 'no' & tissue == FALSE){
# If normalize is true then images are normalized, if tissue is false then we find the tissue mask
## Return both
return(list(mimosa_dataframe = mimosa_dataframe, top_voxels = top_voxels, smoothed = smoothed,
coupling_intercepts = coupling_intercepts, coupling_slopes = coupling_slopes,
normalized = normalized, tissue_mask = tissue_mask))
}
}