resolving and merging pr with JM.

avalcarcel9 · Nov 20, 2019 · 88e65b5 · 88e65b5
2 parents 7a46d76 + 8f6ac7d
commit 88e65b5
Show file tree

Hide file tree

Showing 14 changed files with 176 additions and 278 deletions.
diff --git a/.gitignore b/.gitignore
@@ -4,3 +4,4 @@ doc
 .Rhistory
 .RData
 .Ruserdata
+inst/doc
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: rtapas
 Title: TAPAS: A thresholding approach for probability map automatic segmentation
     automatic segmentation
-Version: 0.0.3
+Version: 0.1.0
 Authors@R: person(given = "Alessandra", family = "Valcarcel", email = "alval@pennmedicine.upenn.edu", role = c("aut","cre"))
 Maintainer: Alessandra Valcarcel <alval@pennmedicine.upenn.edu>
 Description: Trains and makes predictions using the TAPAS
@@ -13,7 +13,6 @@ BugReports: https://github.com/avalcarcel9/rtapas/issues
 LazyData: true
 Depends: R (>= 2.10)
 Imports: 
-    aliviateR (>= 0.0.2),
     doParallel (>= 1.0.14),
     dplyr (>= 0.7.5.9000),
     foreach (>= 1.4.4),
@@ -31,6 +30,7 @@ Imports:
     stringr (>= 1.3.0),
     tibble (>= 1.7.0)
 Suggests:
+    aliviateR (>= 0.0.2),
     covr,
     knitr,
     remotes,

diff --git a/NAMESPACE b/NAMESPACE
@@ -6,7 +6,6 @@ export(tapas_data_par)
 export(tapas_predict)
 export(tapas_predict_par)
 export(tapas_train)
-importFrom(aliviateR,dsc)
 importFrom(doParallel,registerDoParallel)
 importFrom(dplyr,arrange)
 importFrom(dplyr,bind_rows)
@@ -41,6 +40,7 @@ importFrom(mgcv,gam)
 importFrom(mgcv,predict.gam)
 importFrom(neurobase,check_mask)
 importFrom(neurobase,check_nifti)
+importFrom(neurobase,fast_dice)
 importFrom(neurobase,niftiarr)
 importFrom(neuroim,connComp3D)
 importFrom(oro.nifti,is.nifti)

diff --git a/R/tapas_data.R b/R/tapas_data.R
@@ -20,11 +20,10 @@
 #' provide an ID.
 #' @param verbose A `logical` argument to print messages. Set to `TRUE` by default.
 #' @export
-#' @importFrom aliviateR dsc
 #' @importFrom dplyr bind_rows
 #' @importFrom magrittr "%>%"
 #' @importFrom methods as
-#' @importFrom neurobase check_nifti check_mask niftiarr
+#' @importFrom neurobase check_nifti check_mask niftiarr fast_dice
 #' @importFrom neuroim connComp3D
 #' @importFrom oro.nifti is.nifti
 #' @importFrom tibble tibble
@@ -87,15 +86,20 @@ tapas_data <- function(thresholds = seq(from = 0, to = 1, by = 0.01),
   pmap = c(pmap[mask == 1])
 
   # Obtain a matrix of 0/1 after threhsolding at each threshold value
-  pred_lesion = base::sapply(thresholds, function(x) {base::ifelse(pmap > x, 1, 0)})
+  # pred_lesion = base::sapply(thresholds, function(x) {base::ifelse(pmap > x, 1, 0)})
+  pred_lesion = base::sapply(thresholds, function(x) {
+    x = as.numeric(pmap > x)
+    x[is.na(x)] = 0
+    x
+  })
 
   # initialize a results tibble
   results = tibble::tibble(threshold = thresholds,
                            dsc = thresholds,
                            volume =  thresholds,
                            subject_id = subject_id)
 
-  if(verbose == TRUE){
+  if (verbose == TRUE){
     base::message('# Obtaining threshold level information.')
   }
 
@@ -126,14 +130,25 @@ tapas_data <- function(thresholds = seq(from = 0, to = 1, by = 0.01),
     # Return temp_lmask to binary 0/1
     temp_lmask[temp_lmask > 0] = 1
 
+    if (requireNamespace("aliviateR", quietly = TRUE)) {
+      dice_value = aliviateR::dsc(gold_standard = gold_standard, comp_method = temp_lmask)
+    } else {
+      dice_value = neurobase::fast_dice(gold_standard,
+                                        temp_lmask,
+                                        verbose = verbose)
+    }
+
     results = tibble::tibble(threshold = thresholds[j],
-                             dsc = aliviateR::dsc(gold_standard = gold_standard, comp_method = temp_lmask),
+                             dsc = dice_value,
                              volume =  sum(temp_lmask),
                              subject_id = subject_id)
 
   }
 
-  results = base::lapply(1:length(thresholds), calc_dv, temp_lmask = mask, subject_id = subject_id) %>%
+  results = base::lapply(
+    1:length(thresholds),
+    calc_dv, temp_lmask = mask,
+    subject_id = subject_id) %>%
     dplyr::bind_rows()
 
   base::return(results)

diff --git a/R/tapas_train.R b/R/tapas_train.R
@@ -10,6 +10,7 @@
 #' threshold estimate resulting in Sørensen's–Dice coefficient (DSC) greater than or equal to the `dsc_cutoff`
 #' will be included in training the TAPAS model.
 #' @param verbose A `logical` argument to print messages. Set to `TRUE` by default.
+#' @param ... additional arguments to pass to \code{\link[mgcv]{gam}}
 #' @export
 #' @importFrom dplyr arrange bind_rows filter group_by inner_join ungroup mutate row_number select slice summarize
 #' @importFrom gtools inv.logit logit
@@ -102,8 +103,7 @@
 #' tapas_model$train_data
 #' }
 
-tapas_train <- function(data, dsc_cutoff = 0.03, verbose = TRUE){
-
+tapas_train <- function(data, dsc_cutoff = 0.03, verbose = TRUE, ...){
   # Check that verbose is TRUE or FALSE
   if(is.logical(verbose) == FALSE){
     base::stop('# ERROR: verbose must be logical TRUE to return comments throughout the function or FALSE to silence comments.')
@@ -243,7 +243,7 @@ tapas_train <- function(data, dsc_cutoff = 0.03, verbose = TRUE){
 
   # Fit the TAPAS model
   tapas_model = mgcv::gam(formula = gtools::logit(threshold) ~ s(volume),
-                          data = data)
+                          data = data, ...)
 
   if(verbose == TRUE){
     base::message('# Calculating lower and upper bound clamps.')

diff --git a/man/tapas_train.Rd b/man/tapas_train.Rd
diff --git a/tests/testthat/setup_data.R b/tests/testthat/setup_data.R
@@ -0,0 +1,58 @@
+library(rtapas)
+
+in_ci <- function() {
+  nzchar(Sys.getenv("CI"))
+}
+
+grid = 0.01
+if (in_ci()) {
+  grid = 0.02
+}
+
+
+# Data is provided in the rtapas package as arrays. Below we will convert them to nifti objects.
+# Create a list of gold standard manual segmentation
+train_gold_standard_masks = list(gs1 = rtapas::gs1,
+                                 gs2 = rtapas::gs2,
+                                 gs3 = rtapas::gs3,
+                                 gs4 = rtapas::gs4,
+                                 gs5 = rtapas::gs5,
+                                 gs6 = rtapas::gs6,
+                                 gs7 = rtapas::gs7,
+                                 gs8 = rtapas::gs8,
+                                 gs9 = rtapas::gs9,
+                                 gs10 = rtapas::gs10)
+# Convert the gold standard masks to nifti objects
+train_gold_standard_masks = lapply(train_gold_standard_masks, oro.nifti::nifti)
+
+# Make a list of the training probability maps
+train_probability_maps = list(pmap1 = rtapas::pmap1,
+                              pmap2 = rtapas::pmap2,
+                              pmap3 = rtapas::pmap3,
+                              pmap4 = rtapas::pmap4,
+                              pmap5 = rtapas::pmap5,
+                              pmap6 = rtapas::pmap6,
+                              pmap7 = rtapas::pmap7,
+                              pmap8 = rtapas::pmap8,
+                              pmap9 = rtapas::pmap9,
+                              pmap10 = rtapas::pmap10)
+
+# Convert the probability maps to nifti objects
+train_probability_maps = lapply(train_probability_maps, oro.nifti::nifti)
+# Make a list of the brain masks
+train_brain_masks = list(brain_mask1 = rtapas::brain_mask,
+                         brain_mask2 = rtapas::brain_mask,
+                         brain_mask3 = rtapas::brain_mask,
+                         brain_mask4 = rtapas::brain_mask,
+                         brain_mask5 = rtapas::brain_mask,
+                         brain_mask6 = rtapas::brain_mask,
+                         brain_mask7 = rtapas::brain_mask,
+                         brain_mask8 = rtapas::brain_mask,
+                         brain_mask9 = rtapas::brain_mask,
+                         brain_mask10 = rtapas::brain_mask)
+
+# Convert the brain masks to nifti objects
+train_brain_masks = lapply(train_brain_masks, oro.nifti::nifti)
+
+# Specify training IDs
+train_ids = paste0('subject_', 1:length(train_gold_standard_masks))
diff --git a/tests/testthat/test_parallel_wrappers.R b/tests/testthat/test_parallel_wrappers.R
@@ -16,29 +16,31 @@ testthat::test_that("1. Test tapas_data and tapas_data_par produce the same outp
   pmap1 = oro.nifti::nifti(pmap1)
   brain_mask = oro.nifti::nifti(brain_mask)
 
-  train_data = rtapas::tapas_data(thresholds = seq(from = 0, to = 1, by = 0.01),
-                                  pmap = pmap1,
-                                  gold_standard = gs1,
-                                  mask = brain_mask,
-                                  k = 0,
-                                  subject_id = 'subject_1',
-                                  verbose = FALSE)
+  train_data = rtapas::tapas_data(
+    thresholds = seq(from = 0, to = 1, by = grid),
+    pmap = pmap1,
+    gold_standard = gs1,
+    mask = brain_mask,
+    k = 0,
+    subject_id = 'subject_1',
+    verbose = FALSE)
 
   # Make data into a list for parallel version
   gs1 = list(gs1)
   pmap1 = list(pmap1)
   brain_mask = list(brain_mask)
 
-  train_data_par = tapas_data_par(cores = 1,
-                                  thresholds = seq(from = 0, to = 1, by = 0.01),
-                                  pmap = pmap1,
-                                  gold_standard = gs1,
-                                  mask = brain_mask,
-                                  k = 0,
-                                  subject_id = list('subject_1'),
-                                  ret = TRUE,
-                                  outfile = NULL,
-                                  verbose = FALSE)
+  train_data_par = tapas_data_par(
+    cores = 1,
+    thresholds = seq(from = 0, to = 1, by = grid),
+    pmap = pmap1,
+    gold_standard = gs1,
+    mask = brain_mask,
+    k = 0,
+    subject_id = list('subject_1'),
+    ret = TRUE,
+    outfile = NULL,
+    verbose = FALSE)
 
 
   # Test for equality
@@ -47,64 +49,18 @@ testthat::test_that("1. Test tapas_data and tapas_data_par produce the same outp
 
 testthat::test_that("2. Test tapas_predict and tapas_predict_par produce the same output.", {
 
-  # Data is provided in the rtapas package as arrays. Below we will convert them to nifti objects.
-  # Create a list of gold standard manual segmentation
-  train_gold_standard_masks = list(gs1 = gs1,
-                                   gs2 = gs2,
-                                   gs3 = gs3,
-                                   gs4 = gs4,
-                                   gs5 = gs5,
-                                   gs6 = gs6,
-                                   gs7 = gs7,
-                                   gs8 = gs8,
-                                   gs9 = gs9,
-                                   gs10 = gs10)
-  # Convert the gold standard masks to nifti objects
-  train_gold_standard_masks = lapply(train_gold_standard_masks, oro.nifti::nifti)
-
-  # Make a list of the training probability maps
-  train_probability_maps = list(pmap1 = pmap1,
-                                pmap2 = pmap2,
-                                pmap3 = pmap3,
-                                pmap4 = pmap4,
-                                pmap5 = pmap5,
-                                pmap6 = pmap6,
-                                pmap7 = pmap7,
-                                pmap8 = pmap8,
-                                pmap9 = pmap9,
-                                pmap10 = pmap10)
-
-  # Convert the probability maps to nifti objects
-  train_probability_maps = lapply(train_probability_maps, oro.nifti::nifti)
-  # Make a list of the brain masks
-  train_brain_masks = list(brain_mask1 = brain_mask,
-                           brain_mask2 = brain_mask,
-                           brain_mask3 = brain_mask,
-                           brain_mask4 = brain_mask,
-                           brain_mask5 = brain_mask,
-                           brain_mask6 = brain_mask,
-                           brain_mask7 = brain_mask,
-                           brain_mask8 = brain_mask,
-                           brain_mask9 = brain_mask,
-                           brain_mask10 = brain_mask)
-
-  # Convert the brain masks to nifti objects
-  train_brain_masks = lapply(train_brain_masks, oro.nifti::nifti)
-
-  # Specify training IDs
-  train_ids = paste0('subject_', 1:length(train_gold_standard_masks))
-
   # Obtain training data
-  data = rtapas::tapas_data_par(cores = 1,
-                                thresholds = seq(from = 0, to = 1, by = 0.01),
-                                pmap = train_probability_maps,
-                                gold_standard = train_gold_standard_masks,
-                                mask = train_brain_masks,
-                                k = 0,
-                                subject_id = train_ids,
-                                ret = TRUE,
-                                outfile = NULL,
-                                verbose = FALSE)
+  data = rtapas::tapas_data_par(
+    cores = 1,
+    thresholds = seq(from = 0, to = 1, by = grid),
+    pmap = train_probability_maps,
+    gold_standard = train_gold_standard_masks,
+    mask = train_brain_masks,
+    k = 0,
+    subject_id = train_ids,
+    ret = TRUE,
+    outfile = NULL,
+    verbose = FALSE)
 
   # We can now implement the train_tapas function using the data from tapas_data_par
   tapas_model = rtapas::tapas_train(data = data,

diff --git a/tests/testthat/test_tapas_data.R b/tests/testthat/test_tapas_data.R
@@ -16,13 +16,14 @@ testthat::test_that("Test tapas_data run on first subject sample data matches or
   pmap1 = oro.nifti::nifti(pmap1)
   brain_mask = oro.nifti::nifti(brain_mask)
 
-  train_data = rtapas::tapas_data(thresholds = seq(from = 0, to = 1, by = 0.01),
-                                  pmap = pmap1,
-                                  gold_standard = gs1,
-                                  mask = brain_mask,
-                                  k = 0,
-                                  subject_id = 'subject_1',
-                                  verbose = FALSE)
+  train_data = rtapas::tapas_data(
+    thresholds = seq(from = 0, to = 1, by = grid),
+    pmap = pmap1,
+    gold_standard = gs1,
+    mask = brain_mask,
+    k = 0,
+    subject_id = 'subject_1',
+    verbose = FALSE)
 
   # The first run always succeeds
   expect_known_output(train_data, tmp, print = TRUE, update = FALSE)

diff --git a/tests/testthat/test_tapas_data_par.R b/tests/testthat/test_tapas_data_par.R
@@ -16,14 +16,15 @@ testthat::test_that("Test tapas_data_par run on first subject sample data matche
   pmap1 = oro.nifti::nifti(pmap1)
   brain_mask = oro.nifti::nifti(brain_mask)
 
-  train_data_par = rtapas::tapas_data_par(cores = 1,
-                                          thresholds = seq(from = 0, to = 1, by = 0.01),
-                                          pmap = list(pmap1),
-                                          gold_standard = list(gs1),
-                                          mask = list(brain_mask),
-                                          k = 0,
-                                          subject_id = list('subject_1'),
-                                          verbose = FALSE)
+  train_data_par = rtapas::tapas_data_par(
+    cores = 1,
+    thresholds = seq(from = 0, to = 1, by = grid),
+    pmap = list(pmap1),
+    gold_standard = list(gs1),
+    mask = list(brain_mask),
+    k = 0,
+    subject_id = list('subject_1'),
+    verbose = FALSE)
 
   # The first run always succeeds
   expect_known_output(train_data_par, tmp, print = TRUE, update = FALSE)