diff --git a/.Rbuildignore b/.Rbuildignore index be638aa..24ad20c 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -2,3 +2,8 @@ ^\.Rproj\.user$ README\.md ^\.travis\.yml$ +data/.Rapp.history +^refit_sublime\.R$ +^README\.Rmd$ +sublime_train_data\.rda$ +^appveyor\.yml$ diff --git a/.gitignore b/.gitignore index 2e6ca1b..00b1fac 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ .Rproj.user .Rhistory .RData -*.nii.gz \ No newline at end of file +*.nii.gz +data/.Rapp.history diff --git a/.travis.yml b/.travis.yml index e492576..c29ff7b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,22 +1,27 @@ # R for travis: see documentation at https://docs.travis-ci.com/user/languages/r - language: R sudo: false cache: packages -r_check_args: --as-cran + +services: + - xvfb before_script: - - "export DISPLAY=:99.0" - - "sh -e /etc/init.d/xvfb start" - - sleep 3 + - fname=travis_helpers.sh + - wget -O ${fname} http://bit.ly/travis_helpers + - cat ${fname}; source ${fname}; rm ${fname} + - codename=`linux_codename`; + - if [ "${codename}" = "trusty" ]; + then + start_xvfb ; + fi + +os: + - linux + - osx bioc_required: yes use_bioc: yes r: - release - devel -addons: - apt: - packages: - - tcl8.5-dev - - tk8.5-dev \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index fd65482..794cb89 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,8 +2,8 @@ Package: sublime Type: Package Title: Automatic Lesion Incidence Estimation and Detection using Multi-Modality Longitudinal Magnetic Resonance Images -Version: 1.3 -Date: 2016-09-29 +Version: 1.6.3 +Date: 2019-10-22 Authors@R: c(person(given = "Elizabeth M.", "Sweeney", email = "elizabethmargaretsweeney@gmail.com", role = c("aut", "cre")), person(given = "John", "Muschelli", @@ -16,11 +16,16 @@ Description: Creates probability maps of incident and enlarging lesion voxels License: GPL LazyData: true Depends: - oro.nifti + oro.nifti, + R (>= 2.10) Imports: downloader, graphics, stats, grDevices, AnalyzeFMRI -RoxygenNote: 5.0.1.9000 +RoxygenNote: 7.0.2 +Encoding: UTF-8 +Suggests: + testthat (>= 2.1.0), + neurobase diff --git a/NAMESPACE b/NAMESPACE index 401d436..53b5bd8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,10 +3,11 @@ export(SuBLIME_prediction) export(download_data) export(normalize) +export(sublime_prediction) export(voxel_select) -import(downloader) import(oro.nifti) importFrom(AnalyzeFMRI,GaussSmoothArray) +importFrom(downloader,download) importFrom(grDevices,dev.off) importFrom(grDevices,gray) importFrom(grDevices,pdf) diff --git a/R/SuBLIME_model.R b/R/SuBLIME_model.R index 0fc8add..96897cf 100644 --- a/R/SuBLIME_model.R +++ b/R/SuBLIME_model.R @@ -1,9 +1,30 @@ -#' @title Sublime Predictive model +#' @title SuBLIME Predictive model #' #' @description Predictive model for SuBLIME algorithm -#' @usage SuBLIME_model -#' @format An \code{lm} object, but with data an other things, notably +#' @usage sublime_model +#' @format An \code{glm} object, but with data an other things, notably #' qr removed #' @references Sweeney, E. M., et al. "Automatic lesion incidence estimation and detection in multiple sclerosis using multisequence longitudinal MRI." American Journal of Neuroradiology 34.1 (2013): 68-73. -"SuBLIME_model" +"sublime_model" + + +#' @title SuBLIME Predictive model without PD +#' +#' @description Predictive model for SuBLIME algorithm without PD modality +#' NOTE: this may perform much worse than the original SuBLIME model +#' +#' @usage nopd_sublime_model +#' @format An \code{glm} object, but with data an other things, notably +#' qr removed +"nopd_sublime_model" + + +#' @title SuBLIME Predictive model with only T1 and FLAIR +#' +#' @description Predictive model for SuBLIME algorithm with only T1 and FLAIR. +#' NOTE: this may perform much worse than the original SuBLIME model +#' @usage flairt1_sublime_model +#' @format An \code{glm} object, but with data an other things, notably +#' qr removed +"flairt1_sublime_model" diff --git a/R/SuBLIME_prediction.R b/R/SuBLIME_prediction.R index 8778da2..ed4ba8f 100644 --- a/R/SuBLIME_prediction.R +++ b/R/SuBLIME_prediction.R @@ -20,12 +20,17 @@ #' nifti #' @param time_diff Difference in time (in days) between baseline and #' followup, numeric -#' @param baseline_nawm_mask Baseline Normal Appearing white matter mask, either array or class nifti. -#' Will be coerced to logical usign baseline_nawm_mask $> 0$. If NULL, no NAWM normalization +#' @param baseline_nawm_mask Baseline Normal Appearing white matter mask, +#' either array or class nifti. +#' Will be coerced to logical usign baseline_nawm_mask $> 0$. If NULL, +#' no NAWM normalization #' is done (assumes data is already normalized) -#' @param follow_up_nawm_mask Followup Normal Appearing white matter mask, either array or class nifti. -#' Will be coerced to logical usign follow_up_nawm_mask $> 0$. Defaults to baseline_nawm_mask if -#' not specified. If NULL, no NAWM normalization is done (assumes data is already normalized) +#' @param follow_up_nawm_mask Followup Normal Appearing white matter mask, +#' either array or class nifti. +#' Will be coerced to logical usign follow_up_nawm_mask $> 0$. Defaults to +#' baseline_nawm_mask if +#' not specified. If NULL, no NAWM normalization is done (assumes data is +#' already normalized) #' @param brain_mask Brain mask, either array or class nifti. #' Will be #' coerced to logical usign brain_mask $> 0$. #' @param model Model of class \code{\link{lm}} or set of coefficients. @@ -49,27 +54,32 @@ #' download_data() #' modes = c("FLAIR", "PD", "T2", "VolumetricT1") #' modals = paste0(modes, "norm.nii.gz") -#' base_files = system.file(file.path("01/Baseline", modals), package="SuBLIME") -#' base_imgs = lapply(base_files, readNIfTI, reorient=FALSE) -#' f_files = system.file(file.path("01/FollowUp", modals), package="SuBLIME") +#' base_files = system.file("01", "Baseline", modals, +#' package = "sublime") +#' base_imgs = lapply(base_files, readNIfTI, reorient = FALSE) +#' f_files = system.file("01", "FollowUp", modals, package="sublime") #' f_imgs = lapply(f_files, readNIfTI, reorient=FALSE) #' names(base_imgs) = names(f_imgs) = modes -#' baseline_nawm_file = system.file("01/Baseline/nawm.nii.gz", package="SuBLIME") +#' baseline_nawm_file = system.file("01", "Baseline", +#' "nawm.nii.gz", package="sublime") #' baseline_nawm_mask = readNIfTI(baseline_nawm_file, reorient=FALSE) #' baseline_nawm_mask = drop(baseline_nawm_mask) -#' follow_up_nawm_file = system.file("01/FollowUp/nawm.nii.gz", package="SuBLIME") +#' follow_up_nawm_file = system.file("01", "FollowUp", +#' "nawm.nii.gz", package="sublime") #' follow_up_nawm_mask = readNIfTI(follow_up_nawm_file, reorient=FALSE) -#' brain_file = system.file("01/duramask.nii.gz", package="SuBLIME") +#' brain_file = system.file("01", "duramask.nii.gz", package="sublime") #' brain_mask = readNIfTI(brain_file, reorient=FALSE) #' brain_mask = drop(brain_mask) -#' +#' on_cran = !identical(Sys.getenv("NOT_CRAN"), "true") +#' if (on_cran) { #' follow_up_nawm_mask = NULL #' baseline_nawm_mask = NULL +#' } #' smooth.using = "GaussSmoothArray" #' verbose = TRUE #' time_diff = 10 #' voxsel = TRUE -#' model = SuBLIME_model +#' model = sublime_model #' #voxsel.sigma = s.sigma =diag(3,3) #' #s.ksize = 3 #' #voxsel.ksize = 5 @@ -88,9 +98,26 @@ #' brain_mask = brain_mask, #' voxsel = voxsel, #' model = model, plot.imgs= TRUE, -#' pdfname = "~/Dropbox/SuBLIME_Web_Test/01/pckg_diagnostc.pdf" +#' pdfname = file.path(tempdir(), "pckg_diagnostc.pdf") #' ) -#' +#' +#' +#' nopd_outimg = SuBLIME_prediction( +#' baseline_flair = base_imgs[["FLAIR"]], +#' follow_up_flair= f_imgs[["FLAIR"]], +#' baseline_pd = NULL, +#' follow_up_pd = NULL, +#' baseline_t2 = base_imgs[["T2"]], +#' follow_up_t2 = f_imgs[["T2"]], +#' baseline_t1 = base_imgs[["VolumetricT1"]], +#' follow_up_t1 = f_imgs[["VolumetricT1"]], +#' time_diff = time_diff, +#' baseline_nawm_mask = baseline_nawm_mask, +#' brain_mask = brain_mask, +#' voxsel = TRUE, +#' model = sublime::nopd_sublime_model, plot.imgs= TRUE, +#' pdfname = file.path(tempdir(), "pckg_diagnostc.pdf") +#' ) #' names(base_imgs) = paste0("baseline_", c("flair", "pd", "t2", "t1")) #' names(f_imgs) = paste0("follow_up_", c("flair", "pd", "t2", "t1")) #' attach(base_imgs) @@ -101,17 +128,18 @@ #' @importFrom graphics mtext par #' @importFrom AnalyzeFMRI GaussSmoothArray #' @import oro.nifti -SuBLIME_prediction <- function(baseline_flair, follow_up_flair, baseline_pd, - follow_up_pd, baseline_t2, follow_up_t2, baseline_t1, - follow_up_t1, time_diff, baseline_nawm_mask = NULL, - follow_up_nawm_mask = baseline_nawm_mask, brain_mask, - model = SuBLIME::SuBLIME_model, - voxsel = TRUE, - smooth.using = c("GaussSmoothArray", "none"), - voxsel.sigma = diag(3,3), voxsel.ksize = 5, - s.sigma = diag(3,3), s.ksize = 3, - plot.imgs = FALSE, - slice = 90, pdfname="diag.pdf", verbose = TRUE){ +SuBLIME_prediction <- function( + baseline_flair, follow_up_flair, baseline_pd = NULL, + follow_up_pd = NULL, baseline_t2, follow_up_t2, baseline_t1, + follow_up_t1, time_diff, baseline_nawm_mask = NULL, + follow_up_nawm_mask = baseline_nawm_mask, brain_mask, + model = sublime::sublime_model, + voxsel = TRUE, + smooth.using = c("GaussSmoothArray", "none"), + voxsel.sigma = diag(3,3), voxsel.ksize = 5, + s.sigma = diag(3,3), s.ksize = 3, + plot.imgs = FALSE, + slice = 90, pdfname = "diag.pdf", verbose = TRUE){ stopifnot(time_diff > 0) ##requires the package AnalyzeFMRI for volume smoothing## @@ -146,18 +174,33 @@ SuBLIME_prediction <- function(baseline_flair, follow_up_flair, baseline_pd, time_diff = array(time_diff, dim=img.dim) f.imgs = list( - follow_up_flair = follow_up_flair, - follow_up_pd = follow_up_pd, - follow_up_t2 = follow_up_t2, - follow_up_t1 = follow_up_t1 + follow_up_flair = follow_up_flair ) + f.imgs$follow_up_pd = follow_up_pd + f.imgs$follow_up_t2 = follow_up_t2 + f.imgs$follow_up_t1 = follow_up_t1 + b.imgs = list( - baseline_flair = baseline_flair, - baseline_pd = baseline_pd, - baseline_t2 = baseline_t2, - baseline_t1 = baseline_t1 + baseline_flair = baseline_flair, + baseline_t2 = baseline_t2, + baseline_t1 = baseline_t1 ) + b.imgs$baseline_pd = baseline_pd + b.imgs$baseline_t2 = baseline_t2 + b.imgs$baseline_t1 = baseline_t1 + + + f_modes = sub("follow_up", "", names(f.imgs)) + b_modes = sub("baseline", "", names(b.imgs)) + sd = c(setdiff(f_modes, b_modes), setdiff(b_modes, f_modes)) + if (length(sd) > 0) { + stop(paste0("We have one imaging modality in baseline or ", + "followup but not the respective", "followup/base", + "data. Must make sure specifying both or both = NULL")) + } + + #### check image dimensions sapply(f.imgs, function(x){ @@ -176,7 +219,7 @@ SuBLIME_prediction <- function(baseline_flair, follow_up_flair, baseline_pd, norm.b.imgs = lapply(b.imgs, function(image){ x = normalize(image = image, mask = baseline_nawm_mask) }) - + norm.f.imgs = lapply(f.imgs, function(image){ x = normalize(image = image, mask = follow_up_nawm_mask) }) @@ -187,7 +230,7 @@ SuBLIME_prediction <- function(baseline_flair, follow_up_flair, baseline_pd, rm(list=c("norm.b.imgs", "norm.f.imgs")) names(norm.imgs) = paste0('normalized_', names(norm.imgs)) - modes = c("flair", "pd", "t2", "t1") + # modes = c("flair", "pd", "t2", "t1") ### cleanup # rm(list=paste0("baseline_", modes)) @@ -197,58 +240,90 @@ SuBLIME_prediction <- function(baseline_flair, follow_up_flair, baseline_pd, #### Difference images - FLAIR_diff = norm.imgs$normalized_follow_up_flair - norm.imgs$normalized_baseline_flair - PD_diff = norm.imgs$normalized_follow_up_pd - norm.imgs$normalized_baseline_pd - T2_diff = norm.imgs$normalized_follow_up_t2 - norm.imgs$normalized_baseline_t2 - T1_diff = norm.imgs$normalized_follow_up_t1 - norm.imgs$normalized_baseline_t1 + norm_diff = function(x1, x2) { + if (is.null(x1) || is.null(x2)) { + return(NULL) + } + x1 - x2 + } + FLAIR_diff = norm_diff(norm.imgs$normalized_follow_up_flair, + norm.imgs$normalized_baseline_flair) + PD_diff = norm_diff(norm.imgs$normalized_follow_up_pd, + norm.imgs$normalized_baseline_pd) + T2_diff = norm_diff(norm.imgs$normalized_follow_up_t2, + norm.imgs$normalized_baseline_t2) + T1_diff = norm_diff(norm.imgs$normalized_follow_up_t1, + norm.imgs$normalized_baseline_t1) + # FLAIR_diff = + # norm.imgs$normalized_follow_up_flair - norm.imgs$normalized_baseline_flair + # PD_diff = norm.imgs$normalized_follow_up_pd - norm.imgs$normalized_baseline_pd + # T2_diff = norm.imgs$normalized_follow_up_t2 - norm.imgs$normalized_baseline_t2 + # T1_diff = norm.imgs$normalized_follow_up_t1 - norm.imgs$normalized_baseline_t1 + # if (plot.imgs){ plotimage = function(img, name){ + if (is.null(img)) { + return(invisible(NULL)) + } oro.nifti::image(img, col = gray((0:32)/32), xaxt = 'n', yaxt = 'n' ) mtext(name, SOUTH<-1, line=-1.5, adj=.95, cex=1, col="white", outer=FALSE) } - pdfmaker = pdfmaker[1] - pdf(pdfname) - par(mfrow = c(2,4)) - par(mar=c(0, 0, 0, 0)) - - plotimage(norm.imgs$normalized_baseline_flair[,,slice], "F.base") - plotimage(norm.imgs$normalized_baseline_pd[,,slice], "PD.base") - plotimage(norm.imgs$normalized_baseline_t2[,,slice], "T2.base") - plotimage(norm.imgs$normalized_baseline_t1[,,slice], "T1.base") - - plotimage(norm.imgs$normalized_follow_up_flair[,,slice], "F.followup") - plotimage(norm.imgs$normalized_follow_up_pd[,,slice], "PD.followup") - plotimage(norm.imgs$normalized_follow_up_t2[,,slice], "T2.followup") - plotimage(norm.imgs$normalized_follow_up_t1[,,slice], "T1.followup") - - - par(mfrow = c(2,2)) - par(mar=c(0, 0, 0, 0)) + # pdfname = pdfname[1] + pdf(pdfname) + par(mfrow = c(2,length(b_modes))) + par(mar=c(0, 0, 0, 0)) - plotimage(FLAIR_diff[,,slice], "F.diff") - plotimage(PD_diff[,,slice], "PD.diff") - plotimage(T2_diff[,,slice], "T2.diff") - plotimage(T1_diff[,,slice], "T1.diff") + plotimage(norm.imgs$normalized_baseline_flair[,,slice], "F.base") + plotimage(norm.imgs$normalized_baseline_pd[,,slice], "PD.base") + plotimage(norm.imgs$normalized_baseline_t2[,,slice], "T2.base") + plotimage(norm.imgs$normalized_baseline_t1[,,slice], "T1.base") + + plotimage(norm.imgs$normalized_follow_up_flair[,,slice], "F.followup") + plotimage(norm.imgs$normalized_follow_up_pd[,,slice], "PD.followup") + plotimage(norm.imgs$normalized_follow_up_t2[,,slice], "T2.followup") + plotimage(norm.imgs$normalized_follow_up_t1[,,slice], "T1.followup") + + + par(mfrow = c(2,2)) + par(mar=c(0, 0, 0, 0)) + + plotimage(FLAIR_diff[,,slice], "F.diff") + plotimage(PD_diff[,,slice], "PD.diff") + plotimage(T2_diff[,,slice], "T2.diff") + plotimage(T1_diff[,,slice], "T1.diff") } ##create dataframe with images for prediction## SuBLIME_data <- data.frame( FLAIR = c(norm.imgs$normalized_follow_up_flair), - PD = c(norm.imgs$normalized_follow_up_pd), - T2 = c(norm.imgs$normalized_follow_up_t2), - T1 = c(norm.imgs$normalized_follow_up_t1), FLAIR_diff = c(FLAIR_diff), - PD_diff = c(PD_diff), - T2_diff = c(T2_diff), - T1_diff = c(T1_diff), time_diff = c(time_diff)) + + SuBLIME_data$PD = c(norm.imgs$normalized_follow_up_pd) + SuBLIME_data$PD_diff = c(PD_diff) + SuBLIME_data$T2 = c(norm.imgs$normalized_follow_up_t2) + SuBLIME_data$T1 = c(norm.imgs$normalized_follow_up_t1) + + SuBLIME_data$T2_diff = c(T2_diff) + SuBLIME_data$T1_diff = c(T1_diff) + + mult_null = function(x1, x2) { + if (is.null(x1) || is.null(x2)) { + return(NULL) + } + x1 * x2 + } SuBLIME_data$"(Intercept)" = 1 - SuBLIME_data$"FLAIR_diff:time_diff" = SuBLIME_data$time_diff * SuBLIME_data$FLAIR_diff - SuBLIME_data$"time_diff:PD_diff" = SuBLIME_data$time_diff * SuBLIME_data$PD_diff - SuBLIME_data$"time_diff:T2_diff" = SuBLIME_data$time_diff * SuBLIME_data$T2_diff - SuBLIME_data$"time_diff:T1_diff" = SuBLIME_data$time_diff * SuBLIME_data$T1_diff + SuBLIME_data$"FLAIR_diff:time_diff" = + mult_null(SuBLIME_data$time_diff, SuBLIME_data$FLAIR_diff) + SuBLIME_data$"time_diff:PD_diff" = + mult_null(SuBLIME_data$time_diff, SuBLIME_data$PD_diff) + SuBLIME_data$"time_diff:T2_diff" = + mult_null(SuBLIME_data$time_diff, SuBLIME_data$T2_diff) + SuBLIME_data$"time_diff:T1_diff" = + mult_null(SuBLIME_data$time_diff, SuBLIME_data$T1_diff) if (verbose){ message("Making Predictions\n") @@ -260,7 +335,7 @@ SuBLIME_prediction <- function(baseline_flair, follow_up_flair, baseline_pd, preds = predict( object = model, newdata = SuBLIME_data, - type= "response", interval="none", se=FALSE) + type = "response", interval="none", se=FALSE) } else if (inherits(model, "matrix")){ rn = rownames(model) cn = colnames(SuBLIME_data) @@ -271,7 +346,7 @@ SuBLIME_prediction <- function(baseline_flair, follow_up_flair, baseline_pd, stopifnot(all(rn %in% cn)) SuBLIME_data = as.matrix(SuBLIME_data[, rn]) preds = SuBLIME_data %*% model - preds = 1/(1+exp(-preds)) + preds = 1/(1 + exp(-preds)) } else if (inherits(model, "numeric")){ rn = names(model) cn = colnames(SuBLIME_data) @@ -284,83 +359,85 @@ SuBLIME_prediction <- function(baseline_flair, follow_up_flair, baseline_pd, preds = SuBLIME_data %*% t(t(model)) preds = 1/(1+exp(-preds)) } - -SuBLIME_predictions <- array(preds, dim = img.dim) - -if (voxsel){ - if (verbose){ - message("Selecting certain voxels\n") + + SuBLIME_predictions <- array(preds, dim = img.dim) + + if (voxsel){ + if (verbose){ + message("Selecting certain voxels\n") + } + ##Create voxel selection mask## + voxel_select_mask <- voxel_select( + normalized_baseline_t2 = norm.imgs$normalized_baseline_t2, + normalized_follow_up_t2 = norm.imgs$normalized_follow_up_t2, + brain_mask = brain_mask, + sigma = voxsel.sigma, ksize = voxsel.ksize) + + if (plot.imgs){ + ##View voxel selection mask + par(mfrow = c(1,1)) + image(voxel_select_mask[,,slice]) + } + + SuBLIME_predictions = SuBLIME_predictions * voxel_select_mask } - ##Create voxel selection mask## - voxel_select_mask <- voxel_select( - normalized_baseline_t2 = norm.imgs$normalized_baseline_t2, - normalized_follow_up_t2 = norm.imgs$normalized_follow_up_t2, - brain_mask = brain_mask, - sigma= voxsel.sigma, ksize = voxsel.ksize) + ##Apply voxel selection mask to SuBLIME predictions## + SuBLIME_predictions_voxel_select <- SuBLIME_predictions if (plot.imgs){ - ##View voxel selection mask + ##View the predictions## par(mfrow = c(1,1)) - image(voxel_select_mask[,,slice]) + image(SuBLIME_predictions_voxel_select[,,slice]) } - SuBLIME_predictions = SuBLIME_predictions * voxel_select_mask -} -##Apply voxel selection mask to SuBLIME predictions## -SuBLIME_predictions_voxel_select <- SuBLIME_predictions - -if (plot.imgs){ - ##View the predictions## - par(mfrow = c(1,1)) - image(SuBLIME_predictions_voxel_select[,,slice]) -} - - -if (verbose){ - message("Smoothing voxel lesion probabilities\n") -} -##Smooth predictions to incorportate spatial information## -if (smooth.using == "GaussSmoothArray"){ - SuBLIME_predictions_voxel_select_smoothed <- AnalyzeFMRI::GaussSmoothArray( - SuBLIME_predictions_voxel_select, - sigma=s.sigma, - ksize=s.ksize, - mask=brain_mask) -} else if (smooth.using == "FSL") { - stop("Not implemented yet") - -} else if (smooth.using == "none"){ - SuBLIME_predictions_voxel_select_smoothed = SuBLIME_predictions_voxel_select -} else { - stop("Smoothing method not implemented") -} - -if (plot.imgs){ - ##View the smoothed predictions## - par(mfrow = c(1,1)) - image(SuBLIME_predictions_voxel_select_smoothed[,,slice]) -} - - -if (inherits(temp.img, "nifti")){ - temp.img@.Data = SuBLIME_predictions_voxel_select_smoothed - cmax = max(temp.img, na.rm=TRUE) - cmax = ifelse(is.finite(cmax), cmax, 0) - cmin = min(temp.img, na.rm=TRUE) - cmin = ifelse(is.finite(cmin), cmin, 0) - temp.img@cal_max = cmax - temp.img@cal_min = cmin - temp.img@scl_slope = 1 - temp.img@scl_inter = 0 - SuBLIME_predictions_voxel_select_smoothed = temp.img -} - -if (plot.imgs){ + + if (verbose) { + message("Smoothing voxel lesion probabilities\n") + } + ##Smooth predictions to incorportate spatial information## + if (smooth.using == "GaussSmoothArray") { + SuBLIME_predictions_voxel_select_smoothed <- + AnalyzeFMRI::GaussSmoothArray( + SuBLIME_predictions_voxel_select, + sigma = s.sigma, + ksize = s.ksize, + mask = brain_mask + ) + } else if (smooth.using == "FSL") { + stop("Not implemented yet") + + } else if (smooth.using == "none"){ + SuBLIME_predictions_voxel_select_smoothed = SuBLIME_predictions_voxel_select + } else { + stop("Smoothing method not implemented") + } + + if (plot.imgs){ + ##View the smoothed predictions## + par(mfrow = c(1,1)) + image(SuBLIME_predictions_voxel_select_smoothed[,,slice]) + } + + + if (inherits(temp.img, "nifti")){ + temp.img@.Data = SuBLIME_predictions_voxel_select_smoothed + cmax = max(temp.img, na.rm=TRUE) + cmax = ifelse(is.finite(cmax), cmax, 0) + cmin = min(temp.img, na.rm=TRUE) + cmin = ifelse(is.finite(cmin), cmin, 0) + temp.img@cal_max = cmax + temp.img@cal_min = cmin + temp.img@scl_slope = 1 + temp.img@scl_inter = 0 + SuBLIME_predictions_voxel_select_smoothed = temp.img + } + + if (plot.imgs){ dev.off() -} - -##Return SuBLIME predictions##` -return(SuBLIME_predictions_voxel_select_smoothed) + } + + ##Return SuBLIME predictions##` + return(SuBLIME_predictions_voxel_select_smoothed) } @@ -373,18 +450,18 @@ return(SuBLIME_predictions_voxel_select_smoothed) #' @param force Force download of file even if it exists #' @export #' @return Indicator if the file was downloaded and unzipped -#' @import downloader +#' @importFrom downloader download #' @importFrom utils unzip download_data = function( - folder = system.file(package="SuBLIME"), + folder = system.file(package="sublime"), force = FALSE - ){ - - url = file.path("https://github.com/muschellij2/SuBLIME_package", +){ + + url = file.path("https://github.com/muschellij2/SuBLIME", "raw/data/01.zip") destfile = file.path(folder, "01.zip") if (!file.exists(destfile) | force){ - download(url, destfile=destfile) + downloader::download(url, destfile = destfile) } check_file =file.path(folder, "01/Baseline/nawm.nii.gz") if (!file.exists(check_file)){ @@ -393,5 +470,8 @@ download_data = function( } file.exists(check_file) } - - + +#' @rdname SuBLIME_prediction +#' @export +sublime_prediction = SuBLIME_prediction + diff --git a/R/normalize.R b/R/normalize.R index fb98d08..373d296 100644 --- a/R/normalize.R +++ b/R/normalize.R @@ -13,14 +13,26 @@ #'} normalize <- function(image, mask = NULL){ - if (is.null(mask)){ + if (is.null(mask)) { return(image) } + if (is.character(image)) { + image = readNIfTI(image, reorient = FALSE) + } + if (is.character(mask)) { + mask = readNIfTI(mask, reorient = FALSE) + } ### Check dimensions stopifnot(all.equal(dim(mask)[1:3], dim(image)[1:3])) ### Need a logical mask stopifnot(inherits(mask[1], "logical")) + + ### could also do this + # umask = unique(c(mask)) + # umask = as.numeric(umask) + # stopifnot(all(umask %in% c(0, 1, NA))) + # mask = mask != 0 #### Get indices from the mask - faster because subset one time ind = which(mask) diff --git a/README.Rmd b/README.Rmd new file mode 100644 index 0000000..1e21682 --- /dev/null +++ b/README.Rmd @@ -0,0 +1,34 @@ +--- +output: github_document +--- + +[![Build Status](https://travis-ci.org/muschellij2/SuBLIME.svg?branch=master)](https://travis-ci.org/muschellij2/SuBLIME) +[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/sublime)](http://cran.rstudio.com/web/packages/sublime/index.html) +[![](http://cranlogs.r-pkg.org/badges/grand-total/sublime)](http://cran.rstudio.com/web/packages/sublime/index.html) + + + + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.path = "man/figures/README-", + out.width = "100%" +) +``` +# sublime + +`sublime` is an R package that implements the model from [Automatic Lesion Incidence Estimation and Detection in Multiple Sclerosis Using Multisequence Longitudinal MRI](https://dx.doi.org/10.3174%2Fajnr.A3172). + + +## Installation + +You can install sublime from github with: + +```{r gh-installation, eval = FALSE} +# install.packages("devtools") +devtools::install_github("muschellij2/sublime") +``` + + diff --git a/README.md b/README.md index 03e3b16..40192ef 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,23 @@ -[![Build Status](https://travis-ci.org/muschellij2/SuBLIME_package.svg?branch=master)](https://travis-ci.org/muschellij2/SuBLIME_package) -[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/SuBLIME)](http://cran.rstudio.com/web/packages/SuBLIME/index.html) -[![](http://cranlogs.r-pkg.org/badges/grand-total/SuBLIME)](http://cran.rstudio.com/web/packages/SuBLIME/index.html) -SuBLIME is an R package that implements the model from [Automatic Lesion Incidence Estimation and Detection in Multiple Sclerosis Using Multisequence Longitudinal MRI](https://dx.doi.org/10.3174%2Fajnr.A3172). \ No newline at end of file +[![Build +Status](https://travis-ci.org/muschellij2/SuBLIME.svg?branch=master)](https://travis-ci.org/muschellij2/SuBLIME) +[![CRAN\_Status\_Badge](http://www.r-pkg.org/badges/version/sublime)](http://cran.rstudio.com/web/packages/sublime/index.html) +[![](http://cranlogs.r-pkg.org/badges/grand-total/sublime)](http://cran.rstudio.com/web/packages/sublime/index.html) + + + +# sublime + +`sublime` is an R package that implements the model from [Automatic +Lesion Incidence Estimation and Detection in Multiple Sclerosis Using +Multisequence Longitudinal +MRI](https://dx.doi.org/10.3174%2Fajnr.A3172). + +## Installation + +You can install sublime from github with: + +``` r +# install.packages("devtools") +devtools::install_github("muschellij2/sublime") +``` diff --git a/SuBLIME.Rproj b/SuBLIME.Rproj index 01fef02..01db50f 100644 --- a/SuBLIME.Rproj +++ b/SuBLIME.Rproj @@ -15,5 +15,5 @@ LaTeX: pdfLaTeX BuildType: Package PackageUseDevtools: Yes PackageInstallArgs: --no-multiarch -PackageCheckArgs: --as-cran -PackageRoxygenize: rd,collate,namespace +PackageCheckArgs: --as-cran --run-dontrun +PackageRoxygenize: rd,collate,namespace,vignette diff --git a/appveyor.yml b/appveyor.yml new file mode 100644 index 0000000..23f5095 --- /dev/null +++ b/appveyor.yml @@ -0,0 +1,94 @@ +# Neuroconductor Template +# DO NOT CHANGE the "init" and "install" sections below + +# Download script file from GitHub +init: + - rm -f travis-tool.sh + - ps: | + $ErrorActionPreference = "Stop" + Invoke-WebRequest http://raw.github.com/krlmlr/r-appveyor/master/scripts/appveyor-tool.ps1 -OutFile "..\appveyor-tool.ps1" + Import-Module '..\appveyor-tool.ps1' + +platform: x64 + +environment: + global: + NOT_CRAN: true + USE_RTOOLS: true + # R_CHECK_INSTALL_ARGS: "--install-args=--build --no-multiarch " + R_CHECK_ARGS: "--no-multiarch --no-manual " + R_REMOTES_NO_ERRORS_FROM_WARNINGS: true + WARNINGS_ARE_ERRORS: 1 + + +install: + - ps: Bootstrap + + +cache: + - C:\RLibrary + +# Adapt as necessary starting from here +build_script: + - travis-tool.sh install_deps + +test_script: + - choco install ffmpeg & exit 0 + - choco install dcmtk & exit 0 + - choco install libreoffice-fresh & exit 0 + - if not exist "C:\rtm" mkdir C:\rtm + - set TMP=C:/rtm + - set TEMP=C:/rtm + - set TMPDIR=C:/rtm + # need to move to higher dir because CMake path req + - if not exist "C:\run" mkdir C:\run + - cp -avrTf . C:\run + # need the tool a level above + - cp ../travis-tool.sh C:\ + - cd C:\run + - travis-tool.sh run_tests + - ls + - ls ../ + # - cp ../*.zip ./ + - cp -ravTf C:/run %APPVEYOR_BUILD_FOLDER% + - cd %APPVEYOR_BUILD_FOLDER% + - for %%i in (*.zip) do echo %%~ni> zip_filename.txt + - SET /P PKG_TARBALL= 0$. If NULL, no NAWM normalization +\item{baseline_nawm_mask}{Baseline Normal Appearing white matter mask, +either array or class nifti. +Will be coerced to logical usign baseline_nawm_mask $> 0$. If NULL, +no NAWM normalization is done (assumes data is already normalized)} -\item{follow_up_nawm_mask}{Followup Normal Appearing white matter mask, either array or class nifti. -Will be coerced to logical usign follow_up_nawm_mask $> 0$. Defaults to baseline_nawm_mask if -not specified. If NULL, no NAWM normalization is done (assumes data is already normalized)} +\item{follow_up_nawm_mask}{Followup Normal Appearing white matter mask, +either array or class nifti. +Will be coerced to logical usign follow_up_nawm_mask $> 0$. Defaults to +baseline_nawm_mask if +not specified. If NULL, no NAWM normalization is done (assumes data is +already normalized)} \item{brain_mask}{Brain mask, either array or class nifti. Will be #' coerced to logical usign brain_mask $> 0$.} @@ -87,27 +137,32 @@ Takes in MRI images from followup and gets predictions download_data() modes = c("FLAIR", "PD", "T2", "VolumetricT1") modals = paste0(modes, "norm.nii.gz") -base_files = system.file(file.path("01/Baseline", modals), package="SuBLIME") -base_imgs = lapply(base_files, readNIfTI, reorient=FALSE) -f_files = system.file(file.path("01/FollowUp", modals), package="SuBLIME") +base_files = system.file("01", "Baseline", modals, +package = "sublime") +base_imgs = lapply(base_files, readNIfTI, reorient = FALSE) +f_files = system.file("01", "FollowUp", modals, package="sublime") f_imgs = lapply(f_files, readNIfTI, reorient=FALSE) names(base_imgs) = names(f_imgs) = modes -baseline_nawm_file = system.file("01/Baseline/nawm.nii.gz", package="SuBLIME") +baseline_nawm_file = system.file("01", "Baseline", +"nawm.nii.gz", package="sublime") baseline_nawm_mask = readNIfTI(baseline_nawm_file, reorient=FALSE) baseline_nawm_mask = drop(baseline_nawm_mask) -follow_up_nawm_file = system.file("01/FollowUp/nawm.nii.gz", package="SuBLIME") +follow_up_nawm_file = system.file("01", "FollowUp", +"nawm.nii.gz", package="sublime") follow_up_nawm_mask = readNIfTI(follow_up_nawm_file, reorient=FALSE) -brain_file = system.file("01/duramask.nii.gz", package="SuBLIME") +brain_file = system.file("01", "duramask.nii.gz", package="sublime") brain_mask = readNIfTI(brain_file, reorient=FALSE) brain_mask = drop(brain_mask) - +on_cran = !identical(Sys.getenv("NOT_CRAN"), "true") +if (on_cran) { follow_up_nawm_mask = NULL baseline_nawm_mask = NULL +} smooth.using = "GaussSmoothArray" verbose = TRUE time_diff = 10 voxsel = TRUE -model = SuBLIME_model +model = sublime_model #voxsel.sigma = s.sigma =diag(3,3) #s.ksize = 3 #voxsel.ksize = 5 @@ -126,9 +181,26 @@ baseline_nawm_mask = baseline_nawm_mask, brain_mask = brain_mask, voxsel = voxsel, model = model, plot.imgs= TRUE, -pdfname = "~/Dropbox/SuBLIME_Web_Test/01/pckg_diagnostc.pdf" +pdfname = file.path(tempdir(), "pckg_diagnostc.pdf") ) - + + +nopd_outimg = SuBLIME_prediction( + baseline_flair = base_imgs[["FLAIR"]], +follow_up_flair= f_imgs[["FLAIR"]], +baseline_pd = NULL, +follow_up_pd = NULL, +baseline_t2 = base_imgs[["T2"]], +follow_up_t2 = f_imgs[["T2"]], +baseline_t1 = base_imgs[["VolumetricT1"]], +follow_up_t1 = f_imgs[["VolumetricT1"]], +time_diff = time_diff, +baseline_nawm_mask = baseline_nawm_mask, +brain_mask = brain_mask, +voxsel = TRUE, +model = sublime::nopd_sublime_model, plot.imgs= TRUE, +pdfname = file.path(tempdir(), "pckg_diagnostc.pdf") +) names(base_imgs) = paste0("baseline_", c("flair", "pd", "t2", "t1")) names(f_imgs) = paste0("follow_up_", c("flair", "pd", "t2", "t1")) attach(base_imgs) @@ -139,4 +211,3 @@ attach(f_imgs) predict } \keyword{Sublime_prediction} - diff --git a/man/download_data.Rd b/man/download_data.Rd index 9de7d3d..4ea7adf 100644 --- a/man/download_data.Rd +++ b/man/download_data.Rd @@ -4,7 +4,7 @@ \alias{download_data} \title{Download SuBLIME data} \usage{ -download_data(folder = system.file(package = "SuBLIME"), force = FALSE) +download_data(folder = system.file(package = "sublime"), force = FALSE) } \arguments{ \item{folder}{Folder to download the data - usually SuBLIME folder, @@ -18,4 +18,3 @@ Indicator if the file was downloaded and unzipped \description{ Download test data for examples } - diff --git a/man/flairt1_sublime_model.Rd b/man/flairt1_sublime_model.Rd new file mode 100644 index 0000000..cf4d28b --- /dev/null +++ b/man/flairt1_sublime_model.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sublime_model.R +\docType{data} +\name{flairt1_sublime_model} +\alias{flairt1_sublime_model} +\title{SuBLIME Predictive model with only T1 and FLAIR} +\format{An \code{glm} object, but with data an other things, notably +qr removed} +\usage{ +flairt1_sublime_model +} +\description{ +Predictive model for SuBLIME algorithm with only T1 and FLAIR. +NOTE: this may perform much worse than the original SuBLIME model +} +\keyword{datasets} diff --git a/man/nopd_sublime_model.Rd b/man/nopd_sublime_model.Rd new file mode 100644 index 0000000..8658b5e --- /dev/null +++ b/man/nopd_sublime_model.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sublime_model.R +\docType{data} +\name{nopd_sublime_model} +\alias{nopd_sublime_model} +\title{SuBLIME Predictive model without PD} +\format{An \code{glm} object, but with data an other things, notably +qr removed} +\usage{ +nopd_sublime_model +} +\description{ +Predictive model for SuBLIME algorithm without PD modality +NOTE: this may perform much worse than the original SuBLIME model +} +\keyword{datasets} diff --git a/man/normalize.Rd b/man/normalize.Rd index fad36f3..14885f6 100644 --- a/man/normalize.Rd +++ b/man/normalize.Rd @@ -25,4 +25,3 @@ the mask is normal appearing white matter } } \keyword{normalize} - diff --git a/man/SuBLIME_model.Rd b/man/sublime_model.Rd similarity index 63% rename from man/SuBLIME_model.Rd rename to man/sublime_model.Rd index 9d08153..17d2fa0 100644 --- a/man/SuBLIME_model.Rd +++ b/man/sublime_model.Rd @@ -1,13 +1,13 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/SuBLIME_model.R +% Please edit documentation in R/sublime_model.R \docType{data} -\name{SuBLIME_model} -\alias{SuBLIME_model} -\title{Sublime Predictive model} -\format{An \code{lm} object, but with data an other things, notably +\name{sublime_model} +\alias{sublime_model} +\title{SuBLIME Predictive model} +\format{An \code{glm} object, but with data an other things, notably qr removed} \usage{ -SuBLIME_model +sublime_model } \description{ Predictive model for SuBLIME algorithm @@ -16,4 +16,3 @@ Predictive model for SuBLIME algorithm Sweeney, E. M., et al. "Automatic lesion incidence estimation and detection in multiple sclerosis using multisequence longitudinal MRI." American Journal of Neuroradiology 34.1 (2013): 68-73. } \keyword{datasets} - diff --git a/man/voxel_select.Rd b/man/voxel_select.Rd index 6a63b1f..335d953 100644 --- a/man/voxel_select.Rd +++ b/man/voxel_select.Rd @@ -4,8 +4,13 @@ \alias{voxel_select} \title{Select Potential Voxels} \usage{ -voxel_select(normalized_baseline_t2, normalized_follow_up_t2, brain_mask, - sigma = diag(3, 3), ksize = 5) +voxel_select( + normalized_baseline_t2, + normalized_follow_up_t2, + brain_mask, + sigma = diag(3, 3), + ksize = 5 +) } \arguments{ \item{normalized_baseline_t2}{Baseline T2 image, array or object class nifti that} @@ -31,4 +36,3 @@ GaussSmoothArray } \keyword{Selection} \keyword{Voxel} - diff --git a/refit_sublime.R b/refit_sublime.R new file mode 100644 index 0000000..13c1cf6 --- /dev/null +++ b/refit_sublime.R @@ -0,0 +1,65 @@ +rm(list=ls()) +if (FALSE) { + setwd("/dcl01/smart/data/structural/msmri/SuBLIME_Data") + load("SuBLIME_model.Rdata") + model = SuBLIME_model + + train_data = SuBLIME_model$data + train_data$GOLD_Radio = SuBLIME_model$y + + +} else { + load("data/sublime_model.rda") + model = sublime_model + load("sublime_train_data.rda") +} + +mod_func = function(formula) { + update( + model, + formula = GOLD_Radio ~ FLAIR + T2 + T1 + + FLAIR_diff * time_diff + T2_diff * time_diff + + T1_diff * time_diff + ) +} + + +keep_mod = function(model){ + model$y = c() + model$model = c() + model$residuals = c() + model$fitted.values = c() + model$effects = c() + model$qr$qr = c() + model$linear.predictors = c() + model$weights = c() + model$prior.weights = c() + model$data = c() + attr(model$terms,".Environment") = c() + attr(model$formula,".Environment") = c() + model +} + +nopd_sublime_model = mod_func( + formula = GOLD_Radio ~ FLAIR + T2 + T1 + + FLAIR_diff * time_diff + T2_diff * time_diff + + T1_diff * time_diff) + +nopd_sublime_model = keep_mod(nopd_sublime_model) +save(nopd_sublime_model, + file = "nopd_sublime_model.rda", + compress = "xz", + compression_level = 9) + +flairt1_sublime_model = mod_func( + formula = GOLD_Radio ~ FLAIR + T1 + + FLAIR_diff * time_diff + + T1_diff * time_diff) + +flairt1_sublime_model = keep_mod(flairt1_sublime_model) +save(flairt1_sublime_model, + file = "flairt1_sublime_model.rda", + compress = "xz", + compression_level = 9) + + diff --git a/sublime_train_data.rda b/sublime_train_data.rda new file mode 100644 index 0000000..02dec0f Binary files /dev/null and b/sublime_train_data.rda differ diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..6991626 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(sublime) + +test_check("sublime") diff --git a/tests/testthat/brain_mask.nii.gz b/tests/testthat/brain_mask.nii.gz new file mode 100644 index 0000000..e787c76 Binary files /dev/null and b/tests/testthat/brain_mask.nii.gz differ diff --git a/tests/testthat/full_prediction.nii.gz b/tests/testthat/full_prediction.nii.gz new file mode 100644 index 0000000..1193405 Binary files /dev/null and b/tests/testthat/full_prediction.nii.gz differ diff --git a/tests/testthat/test-array_smoother.R b/tests/testthat/test-array_smoother.R new file mode 100644 index 0000000..250e0cf --- /dev/null +++ b/tests/testthat/test-array_smoother.R @@ -0,0 +1,49 @@ +testthat::context("Trying smoothing") + +if (requireNamespace("neurobase", quietly = TRUE)) { + library(neurobase) +} else { + readnii = function(...) { + suppressWarnings({ + nim = oro.nifti::readNIfTI(..., reorient = FALSE) + }) + nim = oro.nifti::drop_img_dim(nim) + nim = oro.nifti::as.nifti(nim) + nim + } + fast_readnii = readnii + check_nifti = function(..., fast = FALSE) { + lapply(..., readnii) + } + +} + +set.seed(20200222) +x = readnii("full_prediction.nii.gz") +arr = array(x, dim = dim(x)) +mask = readnii("brain_mask.nii.gz") + +testthat::test_that("Checking Smoothing", { + testthat::expect_equal(sum(mask), 1302045L) + testthat::expect_equal(sum(arr), 1470.89250452405) + + voxsel.ksize = 5; + s.sigma = diag(3, 3); + s.ksize = 5 + + for (i in 1:5) { + + result = AnalyzeFMRI::GaussSmoothArray( + arr, sigma = s.sigma, ksize = s.ksize) + testthat::expect_equal(sum(result), 1470.89250452405) + testthat::expect_equal(mean(result), 0.000203695608124163) + testthat::expect_false(isTRUE(all.equal(result, arr))) + + result = AnalyzeFMRI::GaussSmoothArray( + arr, sigma = s.sigma, ksize = s.ksize, mask = mask) + testthat::expect_equal(sum(result), 1470.62787152521) + testthat::expect_equal(mean(result), 0.000203658960592501) + # testthat::expect_false(isTRUE(all.equal(result, arr))) + + } +}) diff --git a/tests/testthat/test-prediction.R b/tests/testthat/test-prediction.R new file mode 100644 index 0000000..92b6e99 --- /dev/null +++ b/tests/testthat/test-prediction.R @@ -0,0 +1,182 @@ +testthat::context("Running core functions") + +if (requireNamespace("neurobase", quietly = TRUE)) { + library(neurobase) +} else { + readnii = function(...) { + suppressWarnings({ + nim = oro.nifti::readNIfTI(..., reorient = FALSE) + }) + nim = oro.nifti::drop_img_dim(nim) + nim = oro.nifti::as.nifti(nim) + nim + } + fast_readnii = readnii + check_nifti = function(..., fast = FALSE) { + lapply(..., readnii) + } + +} + +dl = download_data(folder = tempdir()) +modes = c("FLAIR", "PD", "T2", "VolumetricT1") +modals = paste0(modes, "norm.nii.gz") +base_files = file.path(tempdir(), "01", "Baseline", modals) + +testthat::expect_true(all(file.exists(base_files))) +base_imgs = check_nifti(base_files, fast = TRUE) + +f_files = file.path(tempdir(), "01", "FollowUp", modals) +testthat::expect_true(all(file.exists(f_files))) +f_imgs = check_nifti(f_files, fast = TRUE) +names(base_imgs) = names(f_imgs) = modes + +baseline_nawm_file = file.path(tempdir(), "01", "Baseline", + "nawm.nii.gz") +baseline_nawm_mask = fast_readnii(baseline_nawm_file) +follow_up_nawm_file = file.path(tempdir(), "01", "FollowUp", + "nawm.nii.gz") +follow_up_nawm_mask = fast_readnii(follow_up_nawm_file) +brain_file = file.path(tempdir(), "01", "duramask.nii.gz") +brain_mask = fast_readnii(brain_file) + +testthat::test_that("Downloading data", { + testthat::expect_equal(sum(baseline_nawm_mask), 399426) + testthat::expect_equal(sum(follow_up_nawm_mask), 403133) + testthat::expect_equal(sum(brain_mask), 1302045) +}) + + + + +# on_cran = !identical(Sys.getenv("NOT_CRAN"), "true") +# if (on_cran) { +# follow_up_nawm_mask = NULL +# baseline_nawm_mask = NULL +# } + + + +testthat::context("Running Predictions") + +testthat::test_that("Prediction without Smoothing", { + + + verbose = TRUE + time_diff = 10 + voxsel = TRUE + model = sublime_model + + outimg = SuBLIME_prediction( + baseline_flair = base_imgs[["FLAIR"]], + follow_up_flair = f_imgs[["FLAIR"]], + baseline_pd = base_imgs[["PD"]], + follow_up_pd = f_imgs[["PD"]], + baseline_t2 = base_imgs[["T2"]], + follow_up_t2 = f_imgs[["T2"]], + baseline_t1 = base_imgs[["VolumetricT1"]], + follow_up_t1 = f_imgs[["VolumetricT1"]], + time_diff = time_diff, + baseline_nawm_mask = baseline_nawm_mask, + brain_mask = brain_mask, + voxsel = voxsel, + smooth.using = "none", + model = model, plot.imgs = TRUE, + pdfname = file.path(tempdir(), "pckg_diagnostc.pdf") + ) + testthat::expect_equal(sum(outimg), 1470.89250384455) + testthat::expect_equal(max(outimg), 0.999999295985094) + testthat::expect_equal(sum(outimg > 0.5), 1343L) + + + + + + nopd_outimg = SuBLIME_prediction( + baseline_flair = base_imgs[["FLAIR"]], + follow_up_flair = f_imgs[["FLAIR"]], + baseline_pd = NULL, + follow_up_pd = NULL, + baseline_t2 = base_imgs[["T2"]], + follow_up_t2 = f_imgs[["T2"]], + baseline_t1 = base_imgs[["VolumetricT1"]], + follow_up_t1 = f_imgs[["VolumetricT1"]], + time_diff = time_diff, + baseline_nawm_mask = baseline_nawm_mask, + brain_mask = brain_mask, + voxsel = TRUE, + smooth.using = "none", + model = sublime::nopd_sublime_model, plot.imgs = TRUE, + pdfname = file.path(tempdir(), "pckg_diagnostc.pdf") + ) + + testthat::expect_equal(sum(nopd_outimg), 1325.37639250814) + testthat::expect_equal(max(nopd_outimg), 0.999997878227761) + testthat::expect_equal(sum(nopd_outimg > 0.5), 1221L) + + +}) + +testthat::context("Running Smoothing Predictions") + +testthat::test_that("Prediction with Smoothing", { + + verbose = TRUE + time_diff = 10 + voxsel = TRUE + model = sublime_model + + for (i in 1:4) { + print(i) + outimg = SuBLIME_prediction( + baseline_flair = base_imgs[["FLAIR"]], + follow_up_flair = f_imgs[["FLAIR"]], + baseline_pd = base_imgs[["PD"]], + follow_up_pd = f_imgs[["PD"]], + baseline_t2 = base_imgs[["T2"]], + follow_up_t2 = f_imgs[["T2"]], + baseline_t1 = base_imgs[["VolumetricT1"]], + follow_up_t1 = f_imgs[["VolumetricT1"]], + time_diff = time_diff, + baseline_nawm_mask = baseline_nawm_mask, + brain_mask = brain_mask, + voxsel = voxsel, + smooth.using = "GaussSmoothArray", + model = model, plot.imgs = TRUE, + pdfname = file.path(tempdir(), "pckg_diagnostc.pdf") + ) + rout = round(outimg, 10) + + testthat::expect_equal(sum(rout), 1470.3759855156) + testthat::expect_equal(max(rout), 0.9999779077) + testthat::expect_equal(sum(rout > 0.5), 1028L) + + testthat::expect_equal(sum(outimg), 1812.54867963279) + testthat::expect_equal(max(outimg), 0.999996324449074) + testthat::expect_equal(sum(outimg > 0.5), 1239L) + + + nopd_outimg = SuBLIME_prediction( + baseline_flair = base_imgs[["FLAIR"]], + follow_up_flair = f_imgs[["FLAIR"]], + baseline_pd = NULL, + follow_up_pd = NULL, + baseline_t2 = base_imgs[["T2"]], + follow_up_t2 = f_imgs[["T2"]], + baseline_t1 = base_imgs[["VolumetricT1"]], + follow_up_t1 = f_imgs[["VolumetricT1"]], + time_diff = time_diff, + baseline_nawm_mask = baseline_nawm_mask, + brain_mask = brain_mask, + voxsel = voxsel, + smooth.using = "GaussSmoothArray", + + model = sublime::nopd_sublime_model, plot.imgs = TRUE, + pdfname = file.path(tempdir(), "pckg_diagnostc.pdf") + ) + + testthat::expect_equal(sum(nopd_outimg), 1652.33846616358) + testthat::expect_equal(max(nopd_outimg), 0.999988313344169) + testthat::expect_equal(sum(nopd_outimg > 0.5), 1154L) + } +})