R/timeserieswindow2matrix.R

In neuroconductor/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' Time-series image to matrix of time windows around user-specified events
#'
#' Extract a matrix from a time-series image after applying a mask where each
#' row is a space-time vector
#'
#'
#' @param timeseriesmatrix Input timeseriesmatrix from timeseries2matrix
#' @param mask Input mask of type 'antsImage' ... the mask will be replicated
#' into a 4D image of length timewindow + zeropadvalue.
#' @param eventlist time indices for the events to extract
#' @param timewindow n-timepoints around each event, forward in time.
#' @param zeropadvalue pads the mask by this amount fwd and bwd in time.
#' @param spacing optional 4d spacing vector that will impact smoothing
#' parameters
#' @return Success -- an R matrix of dimensions
#' ntimewindow*sizeofnonzerovaluesinmask*nevents \cr
#' @author Avants BB
#' @examples
#'
#' img <- makeImage( c(10,10,10,50) , 0 )
#' mask <- as.antsImage( array( 1 , dim(img)[1:3] ) )
#' mat<-timeseries2matrix( img, mask )
#' mat <- timeserieswindow2matrix( mat , mask, c(1, 4, 5 ), 2, 0 )
#' print( dim(mat$eventmatrix) )
#' print( dim(mat$mask4d) )
#'
#' ##### another approach
#' dim4<-c( 20, 30, 10, 100 )
#' i1<-10:15
#' i2<-11:18
#' i3<-4:8
#' i4<-10:90
#' arr3d<-array(data = 0, dim = dim4[1:3] )
#' arr3d[i1,i2,i3]<-1
#' arr<-array(data = 0, dim =dim4 )
#' for ( t in i4 ) arr[,,,t]<-t
#' nois<-which( arr > 0 )
#' noisv<-rnorm( length(nois) )
#' arr[ nois ]<-arr[ nois ]+noisv*0.0
#' msk <- as.antsImage( arr3d )
#' img <- as.antsImage( arr )
#' mat<-timeseries2matrix( img, msk )
#' eanat<-sparseDecom(  mat, msk, sparseness=0.1, z=0.5,nvecs=2, its=5,cthresh=0, mycoption=1)
#' eanat2<-sparseDecom( mat,  sparseness=0.1,z=0.5,nvecs=2, its=5,cthresh=0, mycoption=1)
#' enomask<-eanat2$eigenanatomyimages[1,]
#' emask<-eanat$eigenanatomyimages[1,]
#' print( enomask[31:40] )
#' print(   emask[31:40] )
#'
#' # same thing with event matrices ....
#' ttt<-timeserieswindow2matrix( mat, msk, c(20,40,60,70) , 6, 0 )
#' tte<-ttt$eventmatrix
#' eanat<-sparseDecom(  tte, ttt$mask4d, sparseness=-0.9, z=0.5,nvecs=2, its=5,cthresh=0, mycoption=1)
#' eanat2<-sparseDecom( tte,  sparseness=-0.9, z=0.5,nvecs=2, its=5,cthresh=0, mycoption=1)
#' enomask<-eanat2$eigenanatomyimages[,1]
#' # back to timematrix
#' tmat<-matrix( enomask,  nrow=6 )
#' # back to image
#' eimg<-antsImageClone( msk )
#' eimg[ msk == 1 ]<-tmat[1,]
#' # convert image space to evec space
#' emat<-eanat2$eigenanatomyimages
#' # convert emat to events FIXME this does not currently work
#' # eavent<-timeserieswindow2matrix( data.matrix(emat) , msk, 1 , 6, 0 )
#' # emask<-eavent$eventmatrix[1,]
#' #############################
#'
#' @export timeserieswindow2matrix
timeserieswindow2matrix <- function(timeseriesmatrix, mask, eventlist, timewindow,
                                    zeropadvalue = 0, spacing = NA) {
  if (length(dim(timeseriesmatrix)) != 2) {
    print("Mask should be of dimensionality 3")
    return(NA)
  }
  if (length(dim(mask)) != 3) {
    print("ERROR: Mask should be of dimensionality 3")
    return(NA)
  }
  # first create the new mask
  mask4d <- makeImage(c(dim(mask), (timewindow + zeropadvalue * 2)), 0)
  mask4d <- as.array(mask4d)
  for (i in (zeropadvalue + 1):(zeropadvalue + timewindow)) {
    maskvec <- mask[mask >= 0]
    mask4d[, , , i] <- maskvec
  }
  mask4d <- as.antsImage(mask4d)
  if (!is.na(spacing)) {
    stopifnot(length(spacing) == 4)
    antsSetSpacing(mask4d, spacing)
  }
  nvox4dmask <- sum(mask4d == 1)
  if (ncol(timeseriesmatrix) * timewindow != nvox4dmask) {
    print("ERROR: ncol(timeseriesmatrix)*timewindow !=  nvox4dmask")
    return(mask4d)
  }
  nevents <- length(eventlist)
  if (min(eventlist) < 0) {
    print("ERROR: min(eventlist) < 0 ")
    return(mask4d)
  }
  if (max(eventlist) > nrow(timeseriesmatrix)) {
    print("ERROR: max(eventlist) > nrow(timeseriesmatrix) ")
    return(mask4d)
  }
  outmat <- matrix(rep(0, nvox4dmask * nevents), nrow = nevents)
  ct <- 1
  for (i in eventlist) {
    maxrow <- i + (timewindow - 1)
    if (maxrow > nrow(timeseriesmatrix))
      maxrow <- nrow(timeseriesmatrix)
    locmat <- as.numeric(t(timeseriesmatrix[i:maxrow, ]))
    if (length(locmat) == nvox4dmask)
      outmat[ct, ] <- locmat
    ct <- ct + 1
  }
  return(list(eventmatrix = outmat, mask4d = mask4d))
}