diff --git a/DESCRIPTION b/DESCRIPTION
index 22798dd..d891ce7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: LINDA
 Title: Lesion Identification with Neighborhood Data Analysis
-Version: 0.3.2
+Version: 0.4.0
 Authors@R: 
     c(
     person(given = "Dorian",
@@ -17,11 +17,14 @@ Authors@R:
 Description: A neuroimaging toolkit for automatic segmentation
     of chronic stroke lesions.
 License: Apache License | file LICENSE
+Depends: R (>= 2.10)
 Imports:
     ANTsR (>= 0.4.6),
     ANTsRCore (>= 0.6.0.1),
     randomForest,
-    magrittr
+    magrittr,
+    stats,
+    methods
 Suggests: 
     covr,
     knitr,
diff --git a/NAMESPACE b/NAMESPACE
index 1942ea6..8cd48fe 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,14 +1,14 @@
 # Generated by roxygen2: do not edit by hand
 
-export(LINDA.mrvnrfs.predict_chunks)
-export(LINDA.mrvnrfs_chunks)
+export(asymmetry_mask)
+export(linda_mrvnrfs.predict_chunks)
 export(linda_predict)
 export(n4_skull_strip)
 export(run_prediction)
+import(randomForest)
 importFrom(ANTsR,abpBrainExtraction)
 importFrom(ANTsR,abpN4)
 importFrom(ANTsR,composeTransformsToField)
-importFrom(ANTsR,mrvnrfs)
 importFrom(ANTsR,reflectImage)
 importFrom(ANTsR,splitMask)
 importFrom(ANTsRCore,antsApplyTransforms)
@@ -25,6 +25,10 @@ importFrom(ANTsRCore,lappend)
 importFrom(ANTsRCore,makeImage)
 importFrom(ANTsRCore,matrixToImages)
 importFrom(ANTsRCore,resampleImage)
+importFrom(ANTsRCore,resampleImageToTarget)
 importFrom(ANTsRCore,smoothImage)
 importFrom(ANTsRCore,thresholdImage)
 importFrom(magrittr,"%>%")
+importFrom(methods,formalArgs)
+importFrom(stats,median)
+importFrom(stats,predict)
diff --git a/R/asymmetry_mask.R b/R/asymmetry_mask.R
index 2fe26ce..71e2474 100644
--- a/R/asymmetry_mask.R
+++ b/R/asymmetry_mask.R
@@ -1,8 +1,20 @@
-asymmetry_mask = function(img,
-                          brain_mask,
-                          reflaxis = 0,
-                          sigma = 2,
-                          verbose = TRUE) {
+#' Asymmetry Mask
+#'
+#' @param img an \code{antsImage} of T1 image
+#' @param brain_mask an \code{antsImage} of brain mask
+#' @param reflaxis Reflection axis
+#' @param sigma Smoothing sigma
+#' @param verbose print diagnostic messages
+#'
+#' @return A list of the reflection and the mask
+#' @export
+#'
+asymmetry_mask = function(
+  img,
+  brain_mask,
+  reflaxis = 0,
+  sigma = 2,
+  verbose = TRUE) {
   # compute asymmetry mask
   print_msg("Computing asymmetry mask...", verbose = verbose)
 
diff --git a/R/compose_warping.R b/R/compose_warping.R
deleted file mode 100644
index 6ece753..0000000
--- a/R/compose_warping.R
+++ /dev/null
@@ -1,26 +0,0 @@
-compose_warping = function() {
-  warpmni = system.file("extdata", "pennTemplate",
-                        "templateToCh2_1Warp.nii.gz",
-                        package = "LINDA",
-                        mustWork = FALSE)
-
-  if (!file.exists(warpmni)) {
-    transforms = sapply(
-      0:2,
-      function(x) {
-        system.file(
-          "extdata", "pennTemplate",
-          paste0("templateToCh2_1Warp_000", x, ".nii.gz"),
-          package = "LINDA",
-          mustWork = TRUE)
-      })
-    res_img = composeTransformsToField(
-      image = mni,
-      transforms = transforms)
-    outfile = tempfile(fileext = ".nii.gz")
-    antsImageWrite(res_img, outfile)
-    return(outfile)
-  } else {
-    return(warpmni)
-  }
-}
diff --git a/R/getLesionFeatures.R b/R/getLesionFeatures.R
index ede4e85..6ce96f1 100644
--- a/R/getLesionFeatures.R
+++ b/R/getLesionFeatures.R
@@ -42,7 +42,6 @@ getLesionFeatures = function(
   #   kmask = bmask,
   #   verbose = verbose)$segmentation
   feats[[1]] = (conavg - kmean) %>% iMath('Normalize')
-
   # FEAT 2: gradient magnitude
   feats[[2]] = img %>% iMath('Grad') %>% iMath('Normalize')
 
@@ -60,6 +59,7 @@ getLesionFeatures = function(
     iMath(truncate, 0.01, 0.99) %>%
     iMath('Normalize')
 
+
   # FEAT 4: kmean
   feats[[4]] = antsImageClone(kmean)
 
@@ -77,5 +77,8 @@ getLesionFeatures = function(
   # FEAT 6: t1 itself
   feats[[6]] = antsImageClone(img)   # iMath(img,'Normalize')
 
+  names(feats) = c("d_controls", "d_controls_grad_mag",
+                  "n4_d_controls", "kmean", "ref_diff",
+                  "t1")
   return(feats)
 }
diff --git a/R/linda_predict_function.R b/R/linda_predict_function.R
index 5498afb..b539972 100644
--- a/R/linda_predict_function.R
+++ b/R/linda_predict_function.R
@@ -18,7 +18,7 @@
 #'
 #' @importFrom ANTsRCore iMath antsImageRead antsImageWrite antsRegistration
 #' @importFrom ANTsRCore resampleImage smoothImage thresholdImage antsImageClone
-#' @importFrom ANTsRCore antsApplyTransforms is.antsImage
+#' @importFrom ANTsRCore antsApplyTransforms is.antsImage resampleImageToTarget
 #' @importFrom ANTsR abpN4 abpBrainExtraction reflectImage
 #' @importFrom ANTsR composeTransformsToField splitMask
 #' @importFrom magrittr %>%
@@ -290,24 +290,36 @@ linda_predict = function(
 
   seg = out3$segmentation
 
-  L$segmentation = seg
+  seg_file = file.path(outdir, 'Prediction3_template.nii.gz')
+  L$segmentation = seg_file
 
-  antsImageWrite(seg, file.path(outdir, 'Prediction3_template.nii.gz'))
+  antsImageWrite(seg, seg_file)
 
   print_msg("Saving 3rd final prediction in native space...",
             verbose = verbose)
+  segnative_file = file.path(outdir, 'Prediction3_native.nii.gz')
   segnative = out3$segmentation_native
 
-  antsImageWrite(segnative,
-                 file.path(outdir, 'Prediction3_native.nii.gz'))
+  L$segmentation_native = segnative_file
+  antsImageWrite(segnative, segnative_file)
 
   graded_map = out3$multi_res_seg$probs[[1]][[4]]
 
+  grad_file = file.path(outdir, 'Prediction3_graded_map.nii.gz')
+  antsImageWrite(graded_map, grad_file)
+
   # save graded map
   probles = resampleImage(graded_map,
                           dim(tempbrain),
                           useVoxels = 1,
                           interpType = 0)
+  # probles = resampleImageToTarget(
+  #   image = graded_map,
+  #   target = tempbrain,
+  #   interpType = "nearestNeighbor")
+
+
+
   problesnative = antsApplyTransforms(
     fixed = simg,
     moving = probles,
@@ -315,27 +327,29 @@ linda_predict = function(
     interpolator = 'Linear',
     verbose = verbose > 1
   )
+  probles_file = file.path(outdir,
+                            'Prediction3_probability_template.nii.gz')
+  L$lesion_probability_template = probles_file
   print_msg("Saving probabilistic prediction in template space...",
             verbose = verbose)
-  antsImageWrite(probles,
-                 file.path(outdir,
-                           'Prediction3_probability_template.nii.gz'))
+  antsImageWrite(probles, probles_file)
 
+  problesnative_file = file.path(
+    outdir, 'Prediction3_probability_native.nii.gz')
+  L$lesion_probability_native = problesnative_file
   print_msg("Saving probabilistic prediction in native space...",
             verbose = verbose)
-  antsImageWrite(problesnative,
-                 file.path(outdir,
-                           'Prediction3_probability_native.nii.gz'))
-
+  antsImageWrite(problesnative, problesnative_file)
 
 
   # save in MNI coordinates
   print_msg("Transferring data in MNI (ch2) space...",
             verbose = verbose)
 
-  warppenn = file.path(outdir, 'Reg3_sub_to_template_warp.nii.gz')
-  affpenn = file.path(outdir, 'Reg3_sub_to_template_affine.mat')
-
+  # warppenn = file.path(outdir, 'Reg3_sub_to_template_warp.nii.gz')
+  warppenn = reg_to_temp_warp
+  # affpenn = file.path(outdir, 'Reg3_sub_to_template_affine.mat')
+  affpenn = reg_to_temp_aff
 
   mni = system.file("extdata", "pennTemplate", "ch2.nii.gz",
                     package = "LINDA",
@@ -391,11 +405,15 @@ linda_predict = function(
 
   print_msg("Saving subject in MNI (ch2) space...",
             verbose = verbose)
-  antsImageWrite(submni, file.path(outdir, 'Subject_in_MNI.nii.gz'))
+  t1_template = file.path(outdir, 'Subject_in_MNI.nii.gz')
+  antsImageWrite(submni, t1_template)
+
+  L$t1_template = t1_template
 
   print_msg("Saving lesion in MNI (ch2) space...",
             verbose = verbose)
-  antsImageWrite(lesmni, file.path(outdir, 'Lesion_in_MNI.nii.gz'))
-  return(lesmni)
-
+  lesion_template = file.path(outdir, 'Lesion_in_MNI.nii.gz')
+  L$lesion_template = lesion_template
+  antsImageWrite(lesmni, lesion_template)
+  return(L)
 }
diff --git a/R/mrvnrfs_chunks.R b/R/mrvnrfs_chunks.R
deleted file mode 100644
index ca4db7c..0000000
--- a/R/mrvnrfs_chunks.R
+++ /dev/null
@@ -1,399 +0,0 @@
-#' multi-res voxelwise neighborhood random forest segmentation learning
-#'
-#' Represents multiscale feature images as a neighborhood and uses the features
-#' to build a random forest segmentation model from an image population
-#'
-#' @param y list of training labels. either an image or numeric value
-#' @param x a list of lists where each list contains feature images
-#' @param labelmask a mask for the features (all in the same image space)
-#' the labelmask defines the number of parallel samples that will be used
-#' per subject sample. two labels will double the number of predictors
-#' contributed from each feature image.
-#' @param rad vector of dimensionality d define nhood radius
-#' @param nsamples (per subject to enter training)
-#' @param ntrees (for the random forest model)
-#' @param multiResSchedule an integer vector defining multi-res levels
-#' @param asFactors boolean - treat the y entries as factors
-#' @return list a 4-list with the rf model, training vector, feature matrix
-#' and the random mask
-#' @author Avants BB, Tustison NJ, Pustina D
-#'
-#' @examples
-#'
-#' mask<-makeImage( c(10,10), 0 )
-#' mask[ 3:6, 3:6 ]<-1
-#' mask[ 5, 5:6]<-2
-#' ilist<-list()
-#' lablist<-list()
-#' inds<-1:50
-#' scl<-0.33 # a noise parameter
-#' for ( predtype in c("label","scalar") )
-#' {
-#' for ( i in inds ) {
-#'   img<-antsImageClone(mask)
-#'   imgb<-antsImageClone(mask)
-#'   limg<-antsImageClone(mask)
-#'   if ( predtype == "label") {  # 4 class prediction
-#'     img[ 3:6, 3:6 ]<-rnorm(16)*scl+(i %% 4)+scl*mean(rnorm(1))
-#'     imgb[ 3:6, 3:6 ]<-rnorm(16)*scl+(i %% 4)+scl*mean(rnorm(1))
-#'     limg[ 3:6, 3:6 ]<-(i %% 4)+1  # the label image is constant
-#'     }
-#'     if ( predtype == "scalar") {
-#'       img[ 3:6, 3:6 ]<-rnorm(16,1)*scl*(i)+scl*mean(rnorm(1))
-#'       imgb[ 3:6, 3:6 ]<-rnorm(16,1)*scl*(i)+scl*mean(rnorm(1))
-#'       limg<-i^2.0  # a real outcome
-#'       }
-#'     ilist[[i]]<-list(img,imgb)  # two features
-#'     lablist[[i]]<-limg
-#'   }
-#' rad<-rep( 1, 2 )
-#' mr <- c(1.5,1)
-#' rfm<-mrvnrfs( lablist , ilist, mask, rad=rad, multiResSchedule=mr,
-#'      asFactors = (  predtype == "label" ) )
-#' rfmresult<-mrvnrfs.predict( rfm$rflist,
-#'      ilist, mask, rad=rad, asFactors=(  predtype == "label" ),
-#'      multiResSchedule=mr )
-#' if ( predtype == "scalar" )
-#'   print( cor( unlist(lablist) , rfmresult$seg ) )
-#' } # end predtype loop
-#'
-#'
-#' @export LINDA.mrvnrfs_chunks
-LINDA.mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
-                     ntrees=500, multiResSchedule=c(4,2,1), asFactors=TRUE,
-                     voxchunk=1000) {
-  library(randomForest)
-  # check y type
-  yisimg<-TRUE
-  if ( typeof(y[[1]]) == "integer" | typeof(y[[1]]) == "double") yisimg<-FALSE
-  rflist<-list()
-  rfct<-1
-
-  mrcount=0
-  for ( mr in multiResSchedule )
-  {
-    mrcount=mrcount+1
-
-    if (mr != 1) {
-      subdim<-round( dim( labelmask ) / mr )
-      subdim[ subdim < 2*rad+1 ] <- ( 2*rad+1 )[  subdim < 2*rad+1 ]
-      submask<-resampleImage( labelmask, subdim, useVoxels=1, interpType=1 )
-    } else { submask = labelmask }
-
-    ysub<-y
-    xsub<-x
-    nfeats<-length(xsub[[1]])
-
-    # resample xsub and ysub
-    if (mr != 1) {
-      for ( i in 1:(length(xsub)) )
-      {
-        if ( yisimg )
-          ysub[[i]]<-resampleImage( y[[i]], subdim, useVoxels=1, interpType=1 )
-
-        xsub[[i]][[1]]<-resampleImage( xsub[[i]][[1]], subdim, useVoxels=1, 0 )
-        if ( nfeats > 1 )
-          for ( k in 2:nfeats )
-          {
-            xsub[[i]][[k]]<-resampleImage( xsub[[i]][[k]], subdim,
-                                           useVoxels=1, 0 )
-          }
-      }
-    }
-
-    # add newprobs from previous run, already correct dimension
-    if ( rfct > 1 )
-    {
-      for ( kk in 1:length(xsub) )
-      {
-        p1<-unlist( xsub[[kk]] )
-        p2<-unlist(newprobs[[kk]])
-        temp<-lappend(  p1 ,  p2  )
-        xsub[[kk]]<-temp
-      }
-      rm(newprobs)
-    }
-
-    nfeats<-length(xsub[[1]])  # update nfeats with newprobs
-
-    # build model for this mr
-    sol<-vwnrfs( ysub, xsub, submask, rad, nsamples, ntrees, asFactors )
-
-
-
-    # apply model, get chunk of probs, put chunk in masterprobs
-    if (mrcount < length(multiResSchedule)) {   # not last mr, calculate probs for next cycle
-      chunk.seq = seq(1, sum(submask>0), by=voxchunk)
-      predtype<-'response'
-      if ( asFactors ) predtype<-'prob'
-
-      # set up probs to fill:
-      masterprobs=list()
-      for (tt1 in 1:length(xsub)) {
-        masterprobs[[tt1]] = list()
-        if (asFactors)
-          nprob = length( unique( c( as.numeric( submask ) ) ) )
-        else nprob=1
-        for (tt2 in 1:nprob) {
-          masterprobs[[tt1]][[tt2]] = submask*0
-        }
-      } # end creating masterprobs
-
-
-      for (ch in 1:length(chunk.seq)) {
-
-        # set end of this chunk
-        if (ch < length(chunk.seq)) { chnxt=chunk.seq[ch+1]-1
-        } else { chnxt=sum(submask>0) }
-
-        # create mask for this chunk
-        temp=which(submask>0, arr.ind=T)[chunk.seq[ch]:chnxt]
-        nnz = submask>0 ; nnz[-temp]=F
-        cropmask = submask+0
-        cropmask[nnz==F] = 0
-
-        # start filling fm
-        testmat<-t(getNeighborhoodInMask( cropmask, cropmask,
-                                          rad, spatial.info=F, boundary.condition='image' ))
-        hdsz<-nrow(testmat) # neighborhood size
-        nent<-nfeats*ncol(testmat)*nrow(testmat)*length(xsub)*1.0
-        fm<-matrix( nrow=(nrow(testmat)*length(xsub)) ,
-                    ncol=ncol(testmat)*nfeats  )
-        rm( testmat )
-
-        seqby<-seq.int( 1, hdsz*length(xsub)+1, by=hdsz )
-        for ( i in 1:(length(xsub)) )
-        {
-          m1<-t(getNeighborhoodInMask( xsub[[i]][[1]], cropmask,
-                                       rad, spatial.info=F, boundary.condition='image' ))
-          if ( nfeats > 1 )
-            for ( k in 2:nfeats )
-            {
-              m2<-t(getNeighborhoodInMask( xsub[[i]][[k]], cropmask,
-                                           rad, spatial.info=F, boundary.condition='image' ))
-              m1<-cbind( m1, m2 )
-            }
-          nxt<-seqby[ i + 1 ]-1
-          fm[ seqby[i]:nxt, ]<-m1
-        } # end filling fm, ready for predict
-
-        probsrf<-t( predict( sol$rfm, newdata=fm, type=predtype ) )
-
-
-        for ( i in 1:(length(xsub)) )
-        {
-          nxt<-seqby[ i + 1 ]-1
-          probsx<-list(submask)
-          if ( asFactors )
-            probsx<-matrixToImages(probsrf[,seqby[i]:nxt],  cropmask )
-          else probsx<-list( makeImage( cropmask, probsrf[,seqby[i]:nxt] ) )
-
-          for (tt1 in 1:length(probsx)) {
-               masterprobs[[i]][[tt1]][cropmask>0] = probsx[[tt1]][cropmask>0]
-          }  # end filling masterprobs
-        }  # end putting probs to images
-
-        message(paste(ch,'of',length(chunk.seq)))
-      } # end chunk loop, masterprobs is complete now
-
-
-      # resample masterprobs back to original resolution
-      newprobs=masterprobs
-      if ( ! all( dim(masterprobs[[1]][[1]] ) == dim(labelmask) ) ) {
-        if (mrcount < length(multiResSchedule)) {  # not last mr, resample to next mr
-          nextdim = round( dim( labelmask ) / multiResSchedule[mrcount+1] )
-          nextdim[ nextdim < 2*rad+1 ] <- ( 2*rad+1 )[  nextdim < 2*rad+1 ]
-          for ( tt1 in 1:length(masterprobs) ) for (tt2 in 1:length(masterprobs[[tt1]]))
-          {
-            newprobs[[tt1]][[tt2]]<-resampleImage( masterprobs[[tt1]][[tt2]], nextdim,
-                                                   useVoxels=1, 0 )
-          }
-
-        } else {  # last mr, resample to labelmask
-          for ( tt1 in 1:length(masterprobs) ) for (tt2 in 1:length(masterprobs[[tt1]]))
-          {
-            newprobs[[tt1]][[tt2]]<-resampleImage( masterprobs[[tt1]][[tt2]], dim(labelmask),
-                                           useVoxels=1, 0 )
-          }
-        }
-      }  # end if that resamples newprobs for next level
-      rm(masterprobs) ; rm(probsrf) ; rm(fm)
-    }  # end if not last mr that computes newprobs
-
-
-    rm(xsub) ; rm(ysub)
-    rflist[[rfct]]<-sol$rfm
-    rfct<-rfct+1
-  } # mr loop
-  return( list(rflist=rflist, randmask=sol$randmask ) )
-}
-
-
-
-#' multi-res voxelwise neighborhood random forest segmentation
-#'
-#' Represents multiscale feature images as a neighborhood and uses the features
-#' to apply a random forest segmentation model to a new image
-#'
-#' @param rflist a list of random forest models from mrvnrfs
-#' @param x a list of lists where each list contains feature images
-#' @param labelmask a mask for the features (all in the same image space)
-#' @param rad vector of dimensionality d define nhood radius
-#' @param multiResSchedule an integer vector defining multi-res levels
-#' @param asFactors boolean - treat the y entries as factors
-#' @return list a 4-list with the rf model, training vector, feature matrix
-#' and the random mask
-#' @author Avants BB, Tustison NJ, Pustina D
-#'
-#' @export LINDA.mrvnrfs.predict_chunks
-#'
-#' @importFrom ANTsRCore getNeighborhoodInMask matrixToImages
-#' @importFrom ANTsRCore lappend imageListToMatrix makeImage
-#' @importFrom ANTsR mrvnrfs splitMask
-LINDA.mrvnrfs.predict_chunks <- function( rflist, x, labelmask, rad=NA,
-                             multiResSchedule=c(4,2,1), asFactors=TRUE,
-                             voxchunk=1000) {
-  library(randomForest)
-  rfct<-1
-  for ( mr in multiResSchedule )
-  {
-    subdim<-round( dim( labelmask ) / mr )
-    subdim[ subdim < 2*rad+1 ] <- ( 2*rad+1 )[  subdim < 2*rad+1 ]
-    submask<-resampleImage( labelmask, subdim, useVoxels=1,
-                            interpType=1 )
-    xsub<-x
-    if ( rfct > 1 )
-    {
-      for ( kk in 1:length(xsub) )
-      {
-        temp<-lappend(  unlist( xsub[[kk]] ) ,  unlist(newprobs[[kk]])  )
-        xsub[[kk]]<-temp
-      }
-    }
-    nfeats<-length(xsub[[1]])
-
-
-    # just resample xsub
-    for ( i in 1:(length(xsub)) )
-    {
-      xsub[[i]][[1]]<-resampleImage( xsub[[i]][[1]], subdim, useVoxels=1, 0 )
-      if ( nfeats > 1 )
-        for ( k in 2:nfeats )
-        {
-          xsub[[i]][[k]]<-resampleImage( xsub[[i]][[k]], subdim,
-                                         useVoxels=1, 0 )
-        }
-    }
-
-
-    predtype<-'response'
-    if ( asFactors ) predtype<-'prob'
-
-    # apply model, get probs, feed them to next level
-    chunk.seq = seq(1, sum(submask>0), by=voxchunk)
-
-    # set up probs to fill, get number of labels from model
-    masterprobs=list()
-    for (tt1 in 1:length(xsub)) {
-      masterprobs[[tt1]] = list()
-      if (asFactors) {
-        nprob = length(levels(rflist[[1]]$y)) + 1
-#         nprob = length( unique( c( as.numeric( submask ) ) ) )
-      } else { nprob=1 }
-      for (tt2 in 1:nprob) {
-        masterprobs[[tt1]][[tt2]] = submask*0
-      }
-    } # end creating masterprobs
-
-
-    chunkmask = splitMask(submask, voxchunk = voxchunk)
-
-
-    for (ch in 1:max(chunkmask)) {
-
-      # # set end of this chunk # removed block after ANTsR reconfig in July 2017
-      # if (ch < length(chunk.seq)) { chnxt=chunk.seq[ch+1]-1
-      # } else { chnxt=sum(submask>0) }
-      #
-      # # create mask for this chunk
-      # temp=which(submask>0)[chunk.seq[ch]:chnxt]
-      # nnz = submask>0 ; nnz[-temp]=F
-      # cropmask = submask+0
-      # cropmask[nnz==F] = 0
-
-      cropmask = thresholdImage(chunkmask, ch, ch)
-
-      # start filling fm
-      testmat<-t(getNeighborhoodInMask( cropmask, cropmask,
-                                        rad, spatial.info=F, boundary.condition='image' ))
-      hdsz<-nrow(testmat) # neighborhood size
-      nent<-nfeats*ncol(testmat)*nrow(testmat)*length(xsub)*1.0
-      fm<-matrix( nrow=(nrow(testmat)*length(xsub)) ,
-                  ncol=ncol(testmat)*nfeats  )
-      rm( testmat )
-
-      seqby<-seq.int( 1, hdsz*length(xsub)+1, by=hdsz )
-      for ( i in 1:(length(xsub)) )
-      {
-        m1<-t(getNeighborhoodInMask( xsub[[i]][[1]], cropmask,
-                                     rad, spatial.info=F, boundary.condition='image' ))
-        if ( nfeats > 1 )
-          for ( k in 2:nfeats )
-          {
-            m2<-t(getNeighborhoodInMask( xsub[[i]][[k]], cropmask,
-                                         rad, spatial.info=F, boundary.condition='image' ))
-            m1<-cbind( m1, m2 )
-          }
-        nxt<-seqby[ i + 1 ]-1
-        fm[ seqby[i]:nxt, ]<-m1
-      } # end filling fm, ready for predict
-
-
-      probs<-t( predict( rflist[[rfct]] ,newdata=fm, type=predtype) )
-
-      for ( i in 1:(length(xsub)) )
-      {
-        nxt<-seqby[ i + 1 ]-1
-        probsx<-list(submask)
-        if ( asFactors ) {
-          probsx<-matrixToImages(probs[,seqby[i]:nxt],  cropmask )
-        } else { probsx<-list( makeImage( cropmask, probs[,seqby[i]:nxt] ) )  }
-
-
-        for (tt1 in 1:length(probsx)) {
-          masterprobs[[i]][[tt1]][cropmask>0] = probsx[[tt1]][cropmask>0]
-        }  # end filling masterprobs
-       }  # end putting probs to images
-
-
-      # resample masterprobs back to original resolution
-      newprobs=masterprobs
-      if ( ! all( dim(masterprobs[[1]][[1]] ) == dim(labelmask) ) )
-      {
-        for ( tt1 in 1:length(masterprobs) ) for (tt2 in 1:length(masterprobs[[tt1]]))
-        {
-          newprobs[[tt1]][[tt2]]<-resampleImage( masterprobs[[tt1]][[tt2]], dim(labelmask),
-                                                 useVoxels=1, 0 )
-        }
-      }
-      message(paste(ch,'of',length(chunk.seq)))
-    }  # end chunk loop
-    rm(masterprobs)
-    rfct<-rfct+1
-  } # mr loop
-
-  # prediction is finished, create segmentation
-  if ( asFactors )
-  {
-    rfseg=list()
-    for (segno in 1:length(newprobs)) {
-      rfseg[[segno]]<-imageListToMatrix( unlist(newprobs[[segno]]) , labelmask )
-      rfseg[[segno]]<-apply( rfseg[[segno]], FUN=which.max, MARGIN=2)
-      rfseg[[segno]]<-makeImage( labelmask , rfseg[[segno]] )
-    }
-    return( list( seg=rfseg, probs=newprobs ) )
-  }
-  rfseg<-apply( imageListToMatrix( unlist(newprobs) ,
-                                   labelmask ), FUN=median, MARGIN=1 )
-  return( list( seg=rfseg, probs=newprobs ) )
-}
diff --git a/R/mrvnrfs_predict_chunks.R b/R/mrvnrfs_predict_chunks.R
new file mode 100644
index 0000000..7466349
--- /dev/null
+++ b/R/mrvnrfs_predict_chunks.R
@@ -0,0 +1,206 @@
+#' multi-res voxelwise neighborhood random forest segmentation
+#'
+#' Represents multiscale feature images as a neighborhood and uses the features
+#' to apply a random forest segmentation model to a new image
+#'
+#' @param rflist a list of random forest models from mrvnrfs
+#' @param x a list of lists where each list contains feature images
+#' @param labelmask a mask for the features (all in the same image space)
+#' @param rad vector of dimensionality d define nhood radius
+#' @param multiResSchedule an integer vector defining multi-res levels
+#' @param asFactors boolean - treat the y entries as factors
+#' @return list a 4-list with the rf model, training vector, feature matrix
+#' and the random mask
+#' @param verbose print diagnostic messages
+#' @param voxchunk value of maximal voxels to predict at once.
+#' This value is used to split the prediction into smaller chunks
+#' such that memory requirements do not become too big
+#' @author Avants BB, Tustison NJ, Pustina D
+#'
+#' @export
+#' @importFrom ANTsRCore lappend makeImage matrixToImages
+#' @importFrom ANTsRCore getNeighborhoodInMask imageListToMatrix
+#' @importFrom ANTsR splitMask
+#' @importFrom stats predict
+#' @import randomForest
+#'
+linda_mrvnrfs.predict_chunks <- function(
+  rflist, x, labelmask, rad=NA,
+  multiResSchedule=c(4,2,1), asFactors=TRUE,
+  voxchunk=1000,
+  verbose = TRUE) {
+  rfct<-1
+  newprobs = NULL
+
+  for ( mr in multiResSchedule )
+  {
+    subdim<-round( dim( labelmask ) / mr )
+    subdim[ subdim < 2*rad+1 ] <- ( 2*rad+1 )[  subdim < 2*rad+1 ]
+    submask<-resampleImage( labelmask, subdim, useVoxels=1,
+                            interpType=1 )
+    xsub<-x
+    if ( rfct > 1 )
+    {
+      for ( kk in 1:length(xsub) )
+      {
+        temp<-lappend(  unlist( xsub[[kk]] ) ,  unlist(newprobs[[kk]])  )
+        xsub[[kk]]<-temp
+      }
+    }
+    nfeats<-length(xsub[[1]])
+
+
+    # just resample xsub
+    for ( i in 1:(length(xsub)) )
+    {
+      xsub[[i]][[1]]<-resampleImage( xsub[[i]][[1]], subdim, useVoxels=1, 0 )
+      # xsub[[i]][[1]] = resampleImageToTarget(
+      #   image = xsub[[i]][[1]],
+      #   target = submask,
+      #   interpType = "linear")
+      if ( nfeats > 1 )
+        for ( k in 2:nfeats )
+        {
+          xsub[[i]][[k]]<-resampleImage( xsub[[i]][[k]], subdim,
+                                         useVoxels=1, 0 )
+          # xsub[[i]][[k]] = resampleImageToTarget(
+          #   image = xsub[[i]][[k]],
+          #   target = submask,
+          #   interpType = "linear")
+        }
+    }
+
+
+    predtype<-'response'
+    if ( asFactors ) predtype<-'prob'
+
+    # apply model, get probs, feed them to next level
+    chunk.seq = seq(1, sum(submask>0), by=voxchunk)
+
+    # set up probs to fill, get number of labels from model
+    masterprobs=list()
+    for (tt1 in 1:length(xsub)) {
+      masterprobs[[tt1]] = list()
+      if (asFactors) {
+        nprob = length(levels(rflist[[1]]$y)) + 1
+        #         nprob = length( unique( c( as.numeric( submask ) ) ) )
+      } else { nprob=1 }
+      for (tt2 in 1:nprob) {
+        masterprobs[[tt1]][[tt2]] = submask*0
+      }
+    } # end creating masterprobs
+
+
+    chunkmask = splitMask(submask, voxchunk = voxchunk)
+
+
+    for (ch in 1:max(chunkmask)) {
+
+      # # set end of this chunk # removed block after ANTsR reconfig in July 2017
+      # if (ch < length(chunk.seq)) { chnxt=chunk.seq[ch+1]-1
+      # } else { chnxt=sum(submask>0) }
+      #
+      # # create mask for this chunk
+      # temp=which(submask>0)[chunk.seq[ch]:chnxt]
+      # nnz = submask>0 ; nnz[-temp]=F
+      # cropmask = submask+0
+      # cropmask[nnz==F] = 0
+
+      cropmask = thresholdImage(chunkmask, ch, ch)
+
+      print_msg("Neighborhood for mask", verbose = verbose)
+
+      # start filling fm
+      testmat<-t(getNeighborhoodInMask( cropmask, cropmask,
+                                        rad, spatial.info=FALSE,
+                                        boundary.condition='image' ))
+      hdsz<-nrow(testmat) # neighborhood size
+      nent<-nfeats*ncol(testmat)*nrow(testmat)*length(xsub)*1.0
+      fm<-matrix( nrow=(nrow(testmat)*length(xsub)) ,
+                  ncol=ncol(testmat)*nfeats  )
+      rm( testmat )
+
+      print_msg("Neighborhood info from features", verbose = verbose)
+
+      seqby<-seq.int( 1, hdsz*length(xsub)+1, by=hdsz )
+      for ( i in 1:(length(xsub)) )
+      {
+        m1<-t(getNeighborhoodInMask(
+          xsub[[i]][[1]], cropmask,
+          rad, spatial.info=FALSE,
+          boundary.condition='image' ))
+        if ( nfeats > 1 )
+          for ( k in 2:nfeats )
+          {
+            m2<-t(getNeighborhoodInMask(
+              xsub[[i]][[k]], cropmask,
+              rad, spatial.info=FALSE,
+              boundary.condition='image' ))
+            m1<-cbind( m1, m2 )
+          }
+        nxt<-seqby[ i + 1 ]-1
+        fm[ seqby[i]:nxt, ]<-m1
+      } # end filling fm, ready for predict
+
+      print_msg("Predicting from RF", verbose = verbose)
+
+      probs<-t( predict( rflist[[rfct]] ,newdata=fm, type=predtype) )
+
+      for ( i in 1:(length(xsub)) )
+      {
+        nxt<-seqby[ i + 1 ]-1
+        probsx<-list(submask)
+        if ( asFactors ) {
+          probsx<-matrixToImages(probs[,seqby[i]:nxt],  cropmask )
+        } else { probsx<-list( makeImage( cropmask, probs[,seqby[i]:nxt] ) )  }
+
+
+        for (tt1 in 1:length(probsx)) {
+          masterprobs[[i]][[tt1]][cropmask>0] = probsx[[tt1]][cropmask>0]
+        }  # end filling masterprobs
+      }  # end putting probs to images
+
+
+      # resample masterprobs back to original resolution
+      newprobs=masterprobs
+      if ( ! all( dim(masterprobs[[1]][[1]] ) == dim(labelmask) ) )
+      {
+        for ( tt1 in 1:length(masterprobs) ) for (tt2 in 1:length(masterprobs[[tt1]]))
+        {
+          newprobs[[tt1]][[tt2]]<-resampleImage(
+            masterprobs[[tt1]][[tt2]], dim(labelmask),
+            useVoxels=1, 0 )
+          # newprobs[[tt1]][[tt2]] = resampleImageToTarget(
+          #   masterprobs[[tt1]][[tt2]],
+          #   target = labelmask,
+          #   interpType = "linear")
+        }
+      }
+      if (verbose) {
+        message(paste(ch,'of',length(chunk.seq)))
+      }
+    }  # end chunk loop
+    rm(masterprobs)
+    rfct<-rfct+1
+    print(rfct)
+  } # mr loop
+
+  # prediction is finished, create segmentation
+  if ( asFactors )
+  {
+    print_msg("Making Segmentations", verbose = verbose)
+
+    rfseg=list()
+    for (segno in 1:length(newprobs)) {
+      rfseg[[segno]]<-imageListToMatrix( unlist(newprobs[[segno]]) , labelmask )
+      rfseg[[segno]]<-apply( rfseg[[segno]], FUN=which.max, MARGIN=2)
+      rfseg[[segno]]<-makeImage( labelmask , rfseg[[segno]] )
+    }
+    return( list( seg=rfseg, probs=newprobs ) )
+  }
+  print_msg("Setting equal to median", verbose = verbose)
+
+  rfseg<-apply( imageListToMatrix( unlist(newprobs) ,
+                                   labelmask ), FUN=median, MARGIN=1 )
+  return( list( seg=rfseg, probs=newprobs ) )
+}
diff --git a/R/n4_skull_strip.R b/R/n4_skull_strip.R
index 24bac7c..7b9c1f0 100644
--- a/R/n4_skull_strip.R
+++ b/R/n4_skull_strip.R
@@ -11,6 +11,8 @@
 #'
 #' @return A list of output images, brain and corrected
 #' @export
+#' @importFrom methods formalArgs
+#' @importFrom stats median
 n4_skull_strip = function(
   file,
   template = system.file("extdata", "pennTemplate", "template.nii.gz",
@@ -66,7 +68,6 @@ n4_skull_strip = function(
     if (!"verbose" %in% formalArgs(abpBrainExtraction)) {
       args$verbose = NULL
     }
-    n4 = do.call(abpN4, args = args)
     bextract = do.call(abpBrainExtraction, args = args)
     rm(submask)
     submask = bextract$bmask * 1
diff --git a/R/run_prediction.R b/R/run_prediction.R
index dbd365b..9060916 100644
--- a/R/run_prediction.R
+++ b/R/run_prediction.R
@@ -33,7 +33,11 @@ run_prediction = function(
   emptyimg = brainmask * 1
   emptyimg = as.array(emptyimg)
   emptyimg[1:91 , , ] = 0
-  template_half_mask = resampleImage(brainmask * emptyimg, voxel_resampling, 0, 1)
+  template_half_mask = resampleImage(
+    brainmask * emptyimg,
+    voxel_resampling,
+    0,
+    1)
   rm(emptyimg)
 
   # fix TruncateIntensity incompatibility with old ANTsR binaries
@@ -75,6 +79,10 @@ run_prediction = function(
       voxel_resampling,
       useVoxels = 0,
       interpType = 0) * template_half_mask
+    # features[[i]] = resampleImageToTarget(
+    #   image = features[[i]],
+    #   target = template_half_mask,
+    #   interpType = "nearestNeighbor") * template_half_mask
   }
 
   # 1st prediction
@@ -89,13 +97,15 @@ run_prediction = function(
   rflist = list(LINDA::rf_model1,
                 LINDA::rf_model2,
                 LINDA::rf_model3)
-  rfm = list(rflist = rflist) # need to keep this for some bad coding in mrvnrfs_chunks
+  # rfm = list(rflist = rflist) # need to keep this for some bad coding in mrvnrfs_chunks
 
 
   mmseg <- suppressMessages(
-    LINDA.mrvnrfs.predict_chunks(
-      rfm$rflist,
+    # ANTsR::mrvnrfs.predict(
+    linda_mrvnrfs.predict_chunks(
+      rflist,
       list(features),
+      # features,
       predlabel.sub,
       rad = rad,
       multiResSchedule = mr,
@@ -109,6 +119,11 @@ run_prediction = function(
                       dim(template_brain),
                       useVoxels = 1,
                       interpType = 1)
+  # seg = resampleImageToTarget(
+  #   image = prediction,
+  #   target = template_brain,
+  #   interpType = "nearestNeighbor")
+
   seg[seg != 4] = 0
   seg[seg == 4] = 1
   segnative = antsApplyTransforms(
diff --git a/man/LINDA.mrvnrfs_chunks.Rd b/man/LINDA.mrvnrfs_chunks.Rd
deleted file mode 100644
index 7405041..0000000
--- a/man/LINDA.mrvnrfs_chunks.Rd
+++ /dev/null
@@ -1,82 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mrvnrfs_chunks.R
-\name{LINDA.mrvnrfs_chunks}
-\alias{LINDA.mrvnrfs_chunks}
-\title{multi-res voxelwise neighborhood random forest segmentation learning}
-\usage{
-LINDA.mrvnrfs_chunks(y, x, labelmask, rad = NA, nsamples = 1,
-  ntrees = 500, multiResSchedule = c(4, 2, 1), asFactors = TRUE,
-  voxchunk = 1000)
-}
-\arguments{
-\item{y}{list of training labels. either an image or numeric value}
-
-\item{x}{a list of lists where each list contains feature images}
-
-\item{labelmask}{a mask for the features (all in the same image space)
-the labelmask defines the number of parallel samples that will be used
-per subject sample. two labels will double the number of predictors
-contributed from each feature image.}
-
-\item{rad}{vector of dimensionality d define nhood radius}
-
-\item{nsamples}{(per subject to enter training)}
-
-\item{ntrees}{(for the random forest model)}
-
-\item{multiResSchedule}{an integer vector defining multi-res levels}
-
-\item{asFactors}{boolean - treat the y entries as factors}
-}
-\value{
-list a 4-list with the rf model, training vector, feature matrix
-and the random mask
-}
-\description{
-Represents multiscale feature images as a neighborhood and uses the features
-to build a random forest segmentation model from an image population
-}
-\examples{
-
-mask<-makeImage( c(10,10), 0 )
-mask[ 3:6, 3:6 ]<-1
-mask[ 5, 5:6]<-2
-ilist<-list()
-lablist<-list()
-inds<-1:50
-scl<-0.33 # a noise parameter
-for ( predtype in c("label","scalar") )
-{
-for ( i in inds ) {
-  img<-antsImageClone(mask)
-  imgb<-antsImageClone(mask)
-  limg<-antsImageClone(mask)
-  if ( predtype == "label") {  # 4 class prediction
-    img[ 3:6, 3:6 ]<-rnorm(16)*scl+(i \%\% 4)+scl*mean(rnorm(1))
-    imgb[ 3:6, 3:6 ]<-rnorm(16)*scl+(i \%\% 4)+scl*mean(rnorm(1))
-    limg[ 3:6, 3:6 ]<-(i \%\% 4)+1  # the label image is constant
-    }
-    if ( predtype == "scalar") {
-      img[ 3:6, 3:6 ]<-rnorm(16,1)*scl*(i)+scl*mean(rnorm(1))
-      imgb[ 3:6, 3:6 ]<-rnorm(16,1)*scl*(i)+scl*mean(rnorm(1))
-      limg<-i^2.0  # a real outcome
-      }
-    ilist[[i]]<-list(img,imgb)  # two features
-    lablist[[i]]<-limg
-  }
-rad<-rep( 1, 2 )
-mr <- c(1.5,1)
-rfm<-mrvnrfs( lablist , ilist, mask, rad=rad, multiResSchedule=mr,
-     asFactors = (  predtype == "label" ) )
-rfmresult<-mrvnrfs.predict( rfm$rflist,
-     ilist, mask, rad=rad, asFactors=(  predtype == "label" ),
-     multiResSchedule=mr )
-if ( predtype == "scalar" )
-  print( cor( unlist(lablist) , rfmresult$seg ) )
-} # end predtype loop
-
-
-}
-\author{
-Avants BB, Tustison NJ, Pustina D
-}
diff --git a/man/asymmetry_mask.Rd b/man/asymmetry_mask.Rd
new file mode 100644
index 0000000..d2083ca
--- /dev/null
+++ b/man/asymmetry_mask.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/asymmetry_mask.R
+\name{asymmetry_mask}
+\alias{asymmetry_mask}
+\title{Asymmetry Mask}
+\usage{
+asymmetry_mask(img, brain_mask, reflaxis = 0, sigma = 2,
+  verbose = TRUE)
+}
+\arguments{
+\item{img}{an \code{antsImage} of T1 image}
+
+\item{brain_mask}{an \code{antsImage} of brain mask}
+
+\item{reflaxis}{Reflection axis}
+
+\item{sigma}{Smoothing sigma}
+
+\item{verbose}{print diagnostic messages}
+}
+\value{
+A list of the reflection and the mask
+}
+\description{
+Asymmetry Mask
+}
diff --git a/man/LINDA.mrvnrfs.predict_chunks.Rd b/man/linda_mrvnrfs.predict_chunks.Rd
similarity index 67%
rename from man/LINDA.mrvnrfs.predict_chunks.Rd
rename to man/linda_mrvnrfs.predict_chunks.Rd
index 855ea66..b23cfe7 100644
--- a/man/LINDA.mrvnrfs.predict_chunks.Rd
+++ b/man/linda_mrvnrfs.predict_chunks.Rd
@@ -1,11 +1,12 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/mrvnrfs_chunks.R
-\name{LINDA.mrvnrfs.predict_chunks}
-\alias{LINDA.mrvnrfs.predict_chunks}
+% Please edit documentation in R/mrvnrfs_predict_chunks.R
+\name{linda_mrvnrfs.predict_chunks}
+\alias{linda_mrvnrfs.predict_chunks}
 \title{multi-res voxelwise neighborhood random forest segmentation}
 \usage{
-LINDA.mrvnrfs.predict_chunks(rflist, x, labelmask, rad = NA,
-  multiResSchedule = c(4, 2, 1), asFactors = TRUE, voxchunk = 1000)
+linda_mrvnrfs.predict_chunks(rflist, x, labelmask, rad = NA,
+  multiResSchedule = c(4, 2, 1), asFactors = TRUE, voxchunk = 1000,
+  verbose = TRUE)
 }
 \arguments{
 \item{rflist}{a list of random forest models from mrvnrfs}
@@ -19,6 +20,12 @@ LINDA.mrvnrfs.predict_chunks(rflist, x, labelmask, rad = NA,
 \item{multiResSchedule}{an integer vector defining multi-res levels}
 
 \item{asFactors}{boolean - treat the y entries as factors}
+
+\item{voxchunk}{value of maximal voxels to predict at once.
+This value is used to split the prediction into smaller chunks
+such that memory requirements do not become too big}
+
+\item{verbose}{print diagnostic messages}
 }
 \value{
 list a 4-list with the rf model, training vector, feature matrix
diff --git a/original/linda_predict.R b/original/linda_predict.R
index 645b182..e658309 100755
--- a/original/linda_predict.R
+++ b/original/linda_predict.R
@@ -9,7 +9,7 @@ cat(paste(format(Sys.time(), "%H:%M") , 'Starting LINDA', ver, "...\n"))
 #' 0.2.3 - fix bug in 'relfaxis'
 #' 0.2.4 - switching mask.lesion1 from graded to binary
 #' 0.2.5 - fixed scriptdir to allow command line call
-#' 0.2.6 - fixed bug in mrvnrfs_chunks.predict after dimfix 
+#' 0.2.6 - fixed bug in mrvnrfs_chunks.predict after dimfix
 #'         from @jeffduda in ANTsR (03/02/2017)
 #'         removed dynamic set of reflaxis, all = 0 now
 #' 0.2.7 - using splitMask for compatibility with new ANTsR
@@ -27,7 +27,7 @@ if (! is.element("randomForest", installed.packages()[,1])) {
   print("Installing missing `randomForest` package")
   install.packages("randomForest")
   library(randomForest)
-} else { 
+} else {
   library(randomForest)
 }
 
@@ -48,7 +48,7 @@ if (! exists('t1')) {
   }
 }
 
-                                                  
+
 # check this is a nifti extension
 if (length(grep(".nii$", t1)) == 0 & length(grep(".nii.gz$", t1)) == 0) {
   stop(
@@ -81,6 +81,7 @@ if ( length( grep('TruncateIntensity' ,iMath(20,'GetOperations'))) != 0 ) {
 if (!exists('scriptdir'))  scriptdir = dirname(sys.frame(1)$ofile)
 source(file.path(scriptdir, 'getLesionFeatures.R'), echo=F)
 source(file.path(scriptdir, 'mrvnrfs_chunks.R'), echo=F)
+source(file.path(scriptdir, 'mrvnrfs_predict_chunks.R'), echo=F)
 
 
 # load template files
@@ -95,7 +96,7 @@ tempmask = antsImageRead(file.path(scriptdir,'pennTemplate','templateBrainMask.n
 cat(paste(format(Sys.time(), "%H:%M") , "Skull stripping... (long process) \n"))
 for (i in 1:2) {
   n4 = abpN4(img = subimg, mask=submask)
-  bextract = abpBrainExtraction(img = subimg, 
+  bextract = abpBrainExtraction(img = subimg,
                                 tem = temp,
                                 temmask = tempmask)
   submask = bextract$bmask*1
@@ -129,22 +130,22 @@ antsImageWrite(mask.lesion1, file.path(lindadir,'Mask.lesion1_asym.nii.gz'))
 resamplevox = c(2,2,2)
 brainmask=iMath(tempmask, 'MD', 1)
 emptyimg = brainmask*1; emptyimg=as.array(emptyimg); emptyimg[ 1:91 , ,]=0
-brainmaskleftbrain = resampleImage(brainmask*emptyimg, resamplevox, 0, 1) 
+brainmaskleftbrain = resampleImage(brainmask*emptyimg, resamplevox, 0, 1)
 rm(emptyimg)
 
 
-################################ 1st registration 
+################################ 1st registration
 cat(paste(format(Sys.time(), "%H:%M") , "Running 1st registration... \n"))
 reg1=antsRegistration(fixed=simg,moving=tempbrain,typeofTransform = 'SyN',mask=mask.lesion1)
-reg1$warpedfixout = reg1$warpedfixout %>% 
-  iMath(truncate,0.01,0.99) %>% 
+reg1$warpedfixout = reg1$warpedfixout %>%
+  iMath(truncate,0.01,0.99) %>%
   iMath('Normalize')
 tempmask=antsApplyTransforms(moving=submask, fixed=tempbrain,transformlist = reg1$invtransforms, interpolator = 'NearestNeighbor')
 
 # prepare features
 cat(paste(format(Sys.time(), "%H:%M") , "Feature calculation... \n"))
 features = getLesionFeatures(reg1$warpedfixout, tempbrain, scriptdir,tempmask,truncate,reflaxis)
-for (i in 1:length(features)) features[[i]] = resampleImage(features[[i]], resamplevox, 
+for (i in 1:length(features)) features[[i]] = resampleImage(features[[i]], resamplevox,
                                                             useVoxels = 0, interpType = 0) * brainmaskleftbrain
 
 # 1st prediction
@@ -164,7 +165,7 @@ cat(paste(format(Sys.time(), "%H:%M") , "Backprojecting 1st prediction... \n"))
 seg = resampleImage(mmseg1$seg[[1]], dim(tempbrain), useVoxels = 1, interpType = 1)
 seg[seg!=4]=0
 seg[seg==4]=1
-segnative = antsApplyTransforms(fixed = simg, moving = seg, 
+segnative = antsApplyTransforms(fixed = simg, moving = seg,
                                 transformlist = reg1$fwdtransforms, interpolator = 'NearestNeighbor')
 mask.lesion2=submask*1
 mask.lesion2[segnative==1]=0
@@ -172,18 +173,18 @@ antsImageWrite(mask.lesion2, file.path(lindadir,'Mask.lesion2.nii.gz'))
 
 
 
-########################### 2nd registration 
+########################### 2nd registration
 cat(paste(format(Sys.time(), "%H:%M") , "Running 2nd registration... \n"))
 reg2=antsRegistration(fixed=simg,moving=tempbrain,typeofTransform = 'SyN',mask=mask.lesion2)
-reg2$warpedfixout = reg2$warpedfixout %>% 
-  iMath(truncate,0.01,0.99) %>% 
+reg2$warpedfixout = reg2$warpedfixout %>%
+  iMath(truncate,0.01,0.99) %>%
   iMath('Normalize')
 tempmask=antsApplyTransforms(moving=submask, fixed=tempbrain,transformlist = reg2$invtransforms, interpolator = 'NearestNeighbor')
 
 # prepare features
 cat(paste(format(Sys.time(), "%H:%M") , "Feature calculation... \n"))
 features = getLesionFeatures(reg2$warpedfixout, tempbrain, scriptdir, tempmask,truncate,reflaxis)
-for (i in 1:length(features)) features[[i]] = resampleImage(features[[i]], resamplevox, 
+for (i in 1:length(features)) features[[i]] = resampleImage(features[[i]], resamplevox,
                                                             useVoxels = 0, interpType = 0) * brainmaskleftbrain
 
 # 2nd prediction
@@ -203,7 +204,7 @@ cat(paste(format(Sys.time(), "%H:%M") , "Backprojecting 2nd prediction... \n"))
 seg = resampleImage(mmseg2$seg[[1]], dim(tempbrain), useVoxels = 1, interpType = 1)
 seg[seg!=4]=0
 seg[seg==4]=1
-segnative = antsApplyTransforms(fixed = simg, moving = seg, 
+segnative = antsApplyTransforms(fixed = simg, moving = seg,
                                 transformlist = reg2$fwdtransforms, interpolator = 'NearestNeighbor')
 mask.lesion3=submask*1
 mask.lesion3[segnative==1]=0
@@ -211,7 +212,7 @@ antsImageWrite(mask.lesion3, file.path(lindadir,'Mask.lesion3.nii.gz'))
 
 
 
-########################### 3rd registration 
+########################### 3rd registration
 cat(paste(format(Sys.time(), "%H:%M") , "Running 3rd registration... (long process)\n"))
 reg3=antsRegistration(fixed=simg,moving=tempbrain,typeofTransform = 'SyNCC',mask=mask.lesion3)
 # save reg3 results
@@ -222,15 +223,15 @@ file.copy(reg3$fwdtransforms[2], file.path(lindadir , 'Reg3_template_to_sub_affi
 file.copy(reg3$invtransforms[1], file.path(lindadir , 'Reg3_sub_to_template_affine.mat'))
 file.copy(reg3$invtransforms[2], file.path(lindadir , 'Reg3_sub_to_template_warp.nii.gz'))
 
-reg3$warpedfixout = reg3$warpedfixout %>% 
-  iMath(truncate,0.01,0.99) %>% 
+reg3$warpedfixout = reg3$warpedfixout %>%
+  iMath(truncate,0.01,0.99) %>%
   iMath('Normalize')
 tempmask=antsApplyTransforms(moving=submask, fixed=tempbrain,transformlist = reg3$invtransforms, interpolator = 'NearestNeighbor')
 
 # prepare features
 cat(paste(format(Sys.time(), "%H:%M") , "Feature calculation... \n"))
 features = getLesionFeatures(reg3$warpedfixout, tempbrain, scriptdir, tempmask,truncate,reflaxis)
-for (i in 1:length(features)) features[[i]] = resampleImage(features[[i]], resamplevox, 
+for (i in 1:length(features)) features[[i]] = resampleImage(features[[i]], resamplevox,
                                                             useVoxels = 0, interpType = 0) * brainmaskleftbrain
 
 # 3rd prediction
@@ -250,7 +251,7 @@ antsImageWrite(mmseg3$seg[[1]], file.path(lindadir,'Prediction3.nii.gz'))
 seg = resampleImage(mmseg3$seg[[1]], dim(tempbrain), useVoxels = 1, interpType = 1)
 seg[seg!=4]=0
 seg[seg==4]=1
-segnative = antsApplyTransforms(fixed = simg, moving = seg, 
+segnative = antsApplyTransforms(fixed = simg, moving = seg,
                                 transformlist = reg3$fwdtransforms, interpolator = 'NearestNeighbor')
 cat(paste(format(Sys.time(), "%H:%M") , "Saving 3rd final prediction in template space... \n"))
 antsImageWrite(seg, file.path(lindadir,'Prediction3_template.nii.gz'))
@@ -259,7 +260,7 @@ antsImageWrite(segnative, file.path(lindadir,'Prediction3_native.nii.gz'))
 
 # save graded map
 probles = resampleImage(mmseg3$probs[[1]][[4]], dim(tempbrain), useVoxels = 1, interpType = 0)
-problesnative = antsApplyTransforms(fixed = simg, moving = probles, 
+problesnative = antsApplyTransforms(fixed = simg, moving = probles,
                                 transformlist = reg3$fwdtransforms, interpolator = 'Linear')
 cat(paste(format(Sys.time(), "%H:%M") , "Saving probabilistic prediction in template space... \n"))
 antsImageWrite(probles, file.path(lindadir,'Prediction3_probability_template.nii.gz'))
@@ -286,4 +287,4 @@ cat(paste(format(Sys.time(), "%H:%M") , "Saving lesion in MNI (ch2) space... \n"
 antsImageWrite(lesmni, file.path(lindadir,'Lesion_in_MNI.nii.gz'))
 
 
-cat('DONE')
\ No newline at end of file
+cat('DONE')
diff --git a/original/mrvnrfs_chunks.R b/original/mrvnrfs_chunks.R
index fd10eb1..b1c38c3 100755
--- a/original/mrvnrfs_chunks.R
+++ b/original/mrvnrfs_chunks.R
@@ -68,29 +68,29 @@ mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
   if ( typeof(y[[1]]) == "integer" | typeof(y[[1]]) == "double") yisimg<-FALSE
   rflist<-list()
   rfct<-1
-  
+
   mrcount=0
   for ( mr in multiResSchedule )
   {
     mrcount=mrcount+1
 
-    if (mr != 1) { 
+    if (mr != 1) {
       subdim<-round( dim( labelmask ) / mr )
       subdim[ subdim < 2*rad+1 ] <- ( 2*rad+1 )[  subdim < 2*rad+1 ]
       submask<-resampleImage( labelmask, subdim, useVoxels=1, interpType=1 )
     } else { submask = labelmask }
-      
+
     ysub<-y
     xsub<-x
     nfeats<-length(xsub[[1]])
-    
+
     # resample xsub and ysub
     if (mr != 1) {
       for ( i in 1:(length(xsub)) )
       {
-        if ( yisimg ) 
+        if ( yisimg )
           ysub[[i]]<-resampleImage( y[[i]], subdim, useVoxels=1, interpType=1 )
-        
+
         xsub[[i]][[1]]<-resampleImage( xsub[[i]][[1]], subdim, useVoxels=1, 0 )
         if ( nfeats > 1 )
           for ( k in 2:nfeats )
@@ -100,7 +100,7 @@ mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
           }
       }
     }
-    
+
     # add newprobs from previous run, already correct dimension
     if ( rfct > 1 )
     {
@@ -113,23 +113,23 @@ mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
       }
       rm(newprobs)
     }
-    
+
     nfeats<-length(xsub[[1]])  # update nfeats with newprobs
-    
+
     # build model for this mr
     sol<-vwnrfs( ysub, xsub, submask, rad, nsamples, ntrees, asFactors )
-    
-    
-    
+
+
+
     # apply model, get chunk of probs, put chunk in masterprobs
     if (mrcount < length(multiResSchedule)) {   # not last mr, calculate probs for next cycle
       chunk.seq = seq(1, sum(submask>0), by=voxchunk)
       predtype<-'response'
       if ( asFactors ) predtype<-'prob'
-      
+
       # set up probs to fill:
       masterprobs=list()
-      for (tt1 in 1:length(xsub)) { 
+      for (tt1 in 1:length(xsub)) {
         masterprobs[[tt1]] = list()
         if (asFactors)
           nprob = length( unique( c( as.numeric( submask ) ) ) )
@@ -138,20 +138,20 @@ mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
           masterprobs[[tt1]][[tt2]] = submask*0
         }
       } # end creating masterprobs
-        
-      
+
+
       for (ch in 1:length(chunk.seq)) {
-        
+
         # set end of this chunk
-        if (ch < length(chunk.seq)) { chnxt=chunk.seq[ch+1]-1 
+        if (ch < length(chunk.seq)) { chnxt=chunk.seq[ch+1]-1
         } else { chnxt=sum(submask>0) }
-        
+
         # create mask for this chunk
         temp=which(submask>0, arr.ind=T)[chunk.seq[ch]:chnxt]
         nnz = submask>0 ; nnz[-temp]=F
         cropmask = submask+0
         cropmask[nnz==F] = 0
-        
+
         # start filling fm
         testmat<-t(getNeighborhoodInMask( cropmask, cropmask,
                                           rad, spatial.info=F, boundary.condition='image' ))
@@ -160,7 +160,7 @@ mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
         fm<-matrix( nrow=(nrow(testmat)*length(xsub)) ,
                     ncol=ncol(testmat)*nfeats  )
         rm( testmat )
-        
+
         seqby<-seq.int( 1, hdsz*length(xsub)+1, by=hdsz )
         for ( i in 1:(length(xsub)) )
         {
@@ -176,27 +176,27 @@ mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
           nxt<-seqby[ i + 1 ]-1
           fm[ seqby[i]:nxt, ]<-m1
         } # end filling fm, ready for predict
-        
+
         probsrf<-t( predict( sol$rfm, newdata=fm, type=predtype ) )
-        
-        
+
+
         for ( i in 1:(length(xsub)) )
         {
           nxt<-seqby[ i + 1 ]-1
           probsx<-list(submask)
-          if ( asFactors ) 
+          if ( asFactors )
             probsx<-matrixToImages(probsrf[,seqby[i]:nxt],  cropmask )
           else probsx<-list( makeImage( cropmask, probsrf[,seqby[i]:nxt] ) )
-          
+
           for (tt1 in 1:length(probsx)) {
                masterprobs[[i]][[tt1]][cropmask>0] = probsx[[tt1]][cropmask>0]
           }  # end filling masterprobs
         }  # end putting probs to images
-        
+
         message(paste(ch,'of',length(chunk.seq)))
       } # end chunk loop, masterprobs is complete now
-  
-      
+
+
       # resample masterprobs back to original resolution
       newprobs=masterprobs
       if ( ! all( dim(masterprobs[[1]][[1]] ) == dim(labelmask) ) ) {
@@ -208,7 +208,7 @@ mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
             newprobs[[tt1]][[tt2]]<-resampleImage( masterprobs[[tt1]][[tt2]], nextdim,
                                                    useVoxels=1, 0 )
           }
-            
+
         } else {  # last mr, resample to labelmask
           for ( tt1 in 1:length(masterprobs) ) for (tt2 in 1:length(masterprobs[[tt1]]))
           {
@@ -219,8 +219,8 @@ mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
       }  # end if that resamples newprobs for next level
       rm(masterprobs) ; rm(probsrf) ; rm(fm)
     }  # end if not last mr that computes newprobs
-    
-    
+
+
     rm(xsub) ; rm(ysub)
     rflist[[rfct]]<-sol$rfm
     rfct<-rfct+1
@@ -230,166 +230,3 @@ mrvnrfs_chunks <- function( y, x, labelmask, rad=NA, nsamples=1,
 
 
 
-#' multi-res voxelwise neighborhood random forest segmentation
-#'
-#' Represents multiscale feature images as a neighborhood and uses the features
-#' to apply a random forest segmentation model to a new image
-#'
-#' @param rflist a list of random forest models from mrvnrfs
-#' @param x a list of lists where each list contains feature images
-#' @param labelmask a mask for the features (all in the same image space)
-#' @param rad vector of dimensionality d define nhood radius
-#' @param multiResSchedule an integer vector defining multi-res levels
-#' @param asFactors boolean - treat the y entries as factors
-#' @return list a 4-list with the rf model, training vector, feature matrix
-#' and the random mask
-#' @author Avants BB, Tustison NJ, Pustina D
-#'
-#' @export mrvnrfs.predict
-mrvnrfs.predict_chunks <- function( rflist, x, labelmask, rad=NA,
-                             multiResSchedule=c(4,2,1), asFactors=TRUE,
-                             voxchunk=1000) {
-  library(randomForest)
-  rfct<-1
-  for ( mr in multiResSchedule )
-  {
-    subdim<-round( dim( labelmask ) / mr )
-    subdim[ subdim < 2*rad+1 ] <- ( 2*rad+1 )[  subdim < 2*rad+1 ]
-    submask<-resampleImage( labelmask, subdim, useVoxels=1,
-                            interpType=1 )
-    xsub<-x
-    if ( rfct > 1 )
-    {
-      for ( kk in 1:length(xsub) )
-      {
-        temp<-lappend(  unlist( xsub[[kk]] ) ,  unlist(newprobs[[kk]])  )
-        xsub[[kk]]<-temp
-      }
-    }
-    nfeats<-length(xsub[[1]])
-    
-    
-    # just resample xsub
-    for ( i in 1:(length(xsub)) )
-    {
-      xsub[[i]][[1]]<-resampleImage( xsub[[i]][[1]], subdim, useVoxels=1, 0 )
-      if ( nfeats > 1 )
-        for ( k in 2:nfeats )
-        {
-          xsub[[i]][[k]]<-resampleImage( xsub[[i]][[k]], subdim,
-                                         useVoxels=1, 0 )
-        }
-    }
-    
-    
-    predtype<-'response'
-    if ( asFactors ) predtype<-'prob'
-    
-    # apply model, get probs, feed them to next level
-    chunk.seq = seq(1, sum(submask>0), by=voxchunk)
-    
-    # set up probs to fill, get number of labels from model
-    masterprobs=list()
-    for (tt1 in 1:length(xsub)) { 
-      masterprobs[[tt1]] = list()
-      if (asFactors) {
-        nprob = length(levels(rfm$rflist[[1]]$y)) + 1
-#         nprob = length( unique( c( as.numeric( submask ) ) ) )
-      } else { nprob=1 }
-      for (tt2 in 1:nprob) {
-        masterprobs[[tt1]][[tt2]] = submask*0
-      }
-    } # end creating masterprobs
-    
-    
-    chunkmask = splitMask(submask, voxchunk = voxchunk)
-    
-    
-    for (ch in 1:max(chunkmask)) {
-      
-      # # set end of this chunk # removed block after ANTsR reconfig in July 2017
-      # if (ch < length(chunk.seq)) { chnxt=chunk.seq[ch+1]-1 
-      # } else { chnxt=sum(submask>0) }
-      # 
-      # # create mask for this chunk
-      # temp=which(submask>0)[chunk.seq[ch]:chnxt]
-      # nnz = submask>0 ; nnz[-temp]=F
-      # cropmask = submask+0
-      # cropmask[nnz==F] = 0
-      
-      cropmask = thresholdImage(chunkmask, ch, ch)
-      
-      # start filling fm
-      testmat<-t(getNeighborhoodInMask( cropmask, cropmask,
-                                        rad, spatial.info=F, boundary.condition='image' ))
-      hdsz<-nrow(testmat) # neighborhood size
-      nent<-nfeats*ncol(testmat)*nrow(testmat)*length(xsub)*1.0
-      fm<-matrix( nrow=(nrow(testmat)*length(xsub)) ,
-                  ncol=ncol(testmat)*nfeats  )
-      rm( testmat )
-      
-      seqby<-seq.int( 1, hdsz*length(xsub)+1, by=hdsz )
-      for ( i in 1:(length(xsub)) )
-      {
-        m1<-t(getNeighborhoodInMask( xsub[[i]][[1]], cropmask,
-                                     rad, spatial.info=F, boundary.condition='image' ))
-        if ( nfeats > 1 )
-          for ( k in 2:nfeats )
-          {
-            m2<-t(getNeighborhoodInMask( xsub[[i]][[k]], cropmask,
-                                         rad, spatial.info=F, boundary.condition='image' ))
-            m1<-cbind( m1, m2 )
-          }
-        nxt<-seqby[ i + 1 ]-1
-        fm[ seqby[i]:nxt, ]<-m1
-      } # end filling fm, ready for predict
-      
-      
-      probs<-t( predict( rflist[[rfct]] ,newdata=fm, type=predtype) )
-      
-      for ( i in 1:(length(xsub)) )
-      {
-        nxt<-seqby[ i + 1 ]-1
-        probsx<-list(submask)
-        if ( asFactors ) {
-          probsx<-matrixToImages(probs[,seqby[i]:nxt],  cropmask )
-        } else { probsx<-list( makeImage( cropmask, probs[,seqby[i]:nxt] ) )  }
-        
-        
-        for (tt1 in 1:length(probsx)) {
-          masterprobs[[i]][[tt1]][cropmask>0] = probsx[[tt1]][cropmask>0]
-        }  # end filling masterprobs
-       }  # end putting probs to images
-      
-      
-      # resample masterprobs back to original resolution
-      newprobs=masterprobs
-      if ( ! all( dim(masterprobs[[1]][[1]] ) == dim(labelmask) ) )
-      {
-        for ( tt1 in 1:length(masterprobs) ) for (tt2 in 1:length(masterprobs[[tt1]]))
-        {
-          newprobs[[tt1]][[tt2]]<-resampleImage( masterprobs[[tt1]][[tt2]], dim(labelmask),
-                                                 useVoxels=1, 0 )
-        }
-      }
-      message(paste(ch,'of',length(chunk.seq)))
-    }  # end chunk loop
-    rm(masterprobs)
-    rfct<-rfct+1
-  } # mr loop
-  
-  # prediction is finished, create segmentation
-  if ( asFactors )
-  {
-    rfseg=list()
-    for (segno in 1:length(newprobs)) {
-      rfseg[[segno]]<-imageListToMatrix( unlist(newprobs[[segno]]) , labelmask )
-      rfseg[[segno]]<-apply( rfseg[[segno]], FUN=which.max, MARGIN=2)
-      rfseg[[segno]]<-makeImage( labelmask , rfseg[[segno]] )
-    }
-    return( list( seg=rfseg, probs=newprobs ) )
-  }
-  rfseg<-apply( imageListToMatrix( unlist(newprobs) ,
-                                   labelmask ), FUN=median, MARGIN=1 )
-  return( list( seg=rfseg, probs=newprobs ) )
-}
\ No newline at end of file
diff --git a/original/mrvnrfs_chunks_predict.R b/original/mrvnrfs_chunks_predict.R
new file mode 100644
index 0000000..30fbbe5
--- /dev/null
+++ b/original/mrvnrfs_chunks_predict.R
@@ -0,0 +1,164 @@
+#' multi-res voxelwise neighborhood random forest segmentation
+#'
+#' Represents multiscale feature images as a neighborhood and uses the features
+#' to apply a random forest segmentation model to a new image
+#'
+#' @param rflist a list of random forest models from mrvnrfs
+#' @param x a list of lists where each list contains feature images
+#' @param labelmask a mask for the features (all in the same image space)
+#' @param rad vector of dimensionality d define nhood radius
+#' @param multiResSchedule an integer vector defining multi-res levels
+#' @param asFactors boolean - treat the y entries as factors
+#' @return list a 4-list with the rf model, training vector, feature matrix
+#' and the random mask
+#' @author Avants BB, Tustison NJ, Pustina D
+#'
+#' @export mrvnrfs.predict
+mrvnrfs.predict_chunks <- function( rflist, x, labelmask, rad=NA,
+                                    multiResSchedule=c(4,2,1), asFactors=TRUE,
+                                    voxchunk=1000) {
+  library(randomForest)
+  rfct<-1
+  for ( mr in multiResSchedule )
+  {
+    subdim<-round( dim( labelmask ) / mr )
+    subdim[ subdim < 2*rad+1 ] <- ( 2*rad+1 )[  subdim < 2*rad+1 ]
+    submask<-resampleImage( labelmask, subdim, useVoxels=1,
+                            interpType=1 )
+    xsub<-x
+    if ( rfct > 1 )
+    {
+      for ( kk in 1:length(xsub) )
+      {
+        temp<-lappend(  unlist( xsub[[kk]] ) ,  unlist(newprobs[[kk]])  )
+        xsub[[kk]]<-temp
+      }
+    }
+    nfeats<-length(xsub[[1]])
+
+
+    # just resample xsub
+    for ( i in 1:(length(xsub)) )
+    {
+      xsub[[i]][[1]]<-resampleImage( xsub[[i]][[1]], subdim, useVoxels=1, 0 )
+      if ( nfeats > 1 )
+        for ( k in 2:nfeats )
+        {
+          xsub[[i]][[k]]<-resampleImage( xsub[[i]][[k]], subdim,
+                                         useVoxels=1, 0 )
+        }
+    }
+
+
+    predtype<-'response'
+    if ( asFactors ) predtype<-'prob'
+
+    # apply model, get probs, feed them to next level
+    chunk.seq = seq(1, sum(submask>0), by=voxchunk)
+
+    # set up probs to fill, get number of labels from model
+    masterprobs=list()
+    for (tt1 in 1:length(xsub)) {
+      masterprobs[[tt1]] = list()
+      if (asFactors) {
+        nprob = length(levels(rflist[[1]]$y)) + 1
+        #         nprob = length( unique( c( as.numeric( submask ) ) ) )
+      } else { nprob=1 }
+      for (tt2 in 1:nprob) {
+        masterprobs[[tt1]][[tt2]] = submask*0
+      }
+    } # end creating masterprobs
+
+
+    chunkmask = splitMask(submask, voxchunk = voxchunk)
+
+
+    for (ch in 1:max(chunkmask)) {
+
+      # # set end of this chunk # removed block after ANTsR reconfig in July 2017
+      # if (ch < length(chunk.seq)) { chnxt=chunk.seq[ch+1]-1
+      # } else { chnxt=sum(submask>0) }
+      #
+      # # create mask for this chunk
+      # temp=which(submask>0)[chunk.seq[ch]:chnxt]
+      # nnz = submask>0 ; nnz[-temp]=F
+      # cropmask = submask+0
+      # cropmask[nnz==F] = 0
+
+      cropmask = thresholdImage(chunkmask, ch, ch)
+
+      # start filling fm
+      testmat<-t(getNeighborhoodInMask( cropmask, cropmask,
+                                        rad, spatial.info=F, boundary.condition='image' ))
+      hdsz<-nrow(testmat) # neighborhood size
+      nent<-nfeats*ncol(testmat)*nrow(testmat)*length(xsub)*1.0
+      fm<-matrix( nrow=(nrow(testmat)*length(xsub)) ,
+                  ncol=ncol(testmat)*nfeats  )
+      rm( testmat )
+
+      seqby<-seq.int( 1, hdsz*length(xsub)+1, by=hdsz )
+      for ( i in 1:(length(xsub)) )
+      {
+        m1<-t(getNeighborhoodInMask( xsub[[i]][[1]], cropmask,
+                                     rad, spatial.info=F, boundary.condition='image' ))
+        if ( nfeats > 1 )
+          for ( k in 2:nfeats )
+          {
+            m2<-t(getNeighborhoodInMask( xsub[[i]][[k]], cropmask,
+                                         rad, spatial.info=F, boundary.condition='image' ))
+            m1<-cbind( m1, m2 )
+          }
+        nxt<-seqby[ i + 1 ]-1
+        fm[ seqby[i]:nxt, ]<-m1
+      } # end filling fm, ready for predict
+
+
+      probs<-t( predict( rflist[[rfct]] ,newdata=fm, type=predtype) )
+
+      for ( i in 1:(length(xsub)) )
+      {
+        nxt<-seqby[ i + 1 ]-1
+        probsx<-list(submask)
+        if ( asFactors ) {
+          probsx<-matrixToImages(probs[,seqby[i]:nxt],  cropmask )
+        } else { probsx<-list( makeImage( cropmask, probs[,seqby[i]:nxt] ) )  }
+
+
+        for (tt1 in 1:length(probsx)) {
+          masterprobs[[i]][[tt1]][cropmask>0] = probsx[[tt1]][cropmask>0]
+        }  # end filling masterprobs
+      }  # end putting probs to images
+
+
+      # resample masterprobs back to original resolution
+      newprobs=masterprobs
+      if ( ! all( dim(masterprobs[[1]][[1]] ) == dim(labelmask) ) )
+      {
+        for ( tt1 in 1:length(masterprobs) ) for (tt2 in 1:length(masterprobs[[tt1]]))
+        {
+          newprobs[[tt1]][[tt2]]<-resampleImage( masterprobs[[tt1]][[tt2]], dim(labelmask),
+                                                 useVoxels=1, 0 )
+        }
+      }
+      message(paste(ch,'of',length(chunk.seq)))
+    }  # end chunk loop
+    rm(masterprobs)
+    rfct<-rfct+1
+    print(rfct)
+  } # mr loop
+
+  # prediction is finished, create segmentation
+  if ( asFactors )
+  {
+    rfseg=list()
+    for (segno in 1:length(newprobs)) {
+      rfseg[[segno]]<-imageListToMatrix( unlist(newprobs[[segno]]) , labelmask )
+      rfseg[[segno]]<-apply( rfseg[[segno]], FUN=which.max, MARGIN=2)
+      rfseg[[segno]]<-makeImage( labelmask , rfseg[[segno]] )
+    }
+    return( list( seg=rfseg, probs=newprobs ) )
+  }
+  rfseg<-apply( imageListToMatrix( unlist(newprobs) ,
+                                   labelmask ), FUN=median, MARGIN=1 )
+  return( list( seg=rfseg, probs=newprobs ) )
+}
diff --git a/tests/testthat.R b/tests/testthat.R
deleted file mode 100644
index f8d0ed8..0000000
--- a/tests/testthat.R
+++ /dev/null
@@ -1,4 +0,0 @@
-library(testthat)
-library(LINDA)
-
-test_check("LINDA")