R/fmri_simulate_func.R

In TCIU: Spacekime Analytics, Time Complexity and Inferential Uncertainty

#' @title real-valued fMRI data simulation
#' @description a real-valued fMRI data simulation function, used to simply generate a 3D fMRI data associated with brain area with activated parts inside. 
#'
#' @param dim_data a vector of length 3 to identify the dimension of fMRI data user wants to simulate
#' @param ons a vector of the start time points of the time period when the fMRI data receives stimulation
#' @param dur a vector of the time period when the fMRI data receives stimulation
#' @param mask a 3D array of 1’s and 0’s or NULL. To specify the area inside the brain shell.
#' One may use the mask data provided by this package, or generate a 3D array of 1’s and 0’s of the same dimension with the fMRI data to be generated.
#' If NULL, then the function would generate a 3D sphere mask.
#'
#' @details
#' The function \code{fmri_simulate_func} is used to simulate fMRI data with specified dimension and total time points. 
#' The fMRI data can be brain-shaped by using the mask data provided in our package, if the dimension fits the same as our data (c(64, 64. 40)).
#' Otherwise, the function will generate a 3D sphere data with multiple activated part inside. The activated parts can be detected based on the p values.
#' 
#'
#' @author SOCR team <\url{http://socr.umich.edu/people/}>
#'
#' @return an array with the specified dimension
#' @return a list of four elements
#' \itemize{
#'   \item fmri_data - the fMRI data generated by the function as specialized values.
#'   \item mask - mask of the fMRI data.
#'   \item ons - a vector of the start time points of the time period when the fMRI data receives stimulation.
#'   \item dur - a vector of the time period when the fMRI data receives stimulation. Notice that the length of ons 
#'   is equal to the length of dur, and all the time period when the data does not receive the simulations have the same 
#'   duration as its former 'on' time period.
#'   \item on_time - a vector that specifies when motion happens.
#' }
#' @export
#'
#' @examples
#' # sample 3D data of mask provided by the package
#' dim(mask)
#' 
#' # the input dimension is the dimension we want for our simulated fMRI data
#' fmri_generate = fmri_simulate_func(dim_data = c(64, 64, 40), mask = mask, 
#'                                    ons = c(1, 21, 41, 61, 81, 101, 121, 141), 
#'                                    dur = c(10, 10, 10, 10, 10, 10, 10, 10))
#' 

fmri_simulate_func = function(dim_data,
                              mask = NULL,
                              ons = c(1, 21, 41, 61, 81, 101, 121, 141),
                              dur = c(10, 10, 10, 10, 10, 10, 10, 10)) {
    
    dimx = dim_data[1]
    dimy = dim_data[2]
    dimz = dim_data[3]
    motor_3d = array(1, dim = dim_data)
    r_motor = (min(dimx, dimy, dimz) + 1)%/%11
    
    
    # mask is provided
    if (is.null(mask) == FALSE) {
        
        for (x in 1:dimx) {
            for (y in 1:dimy) {
                for (z in 1:dimz) {
                  x0 = (1 + dimx)%/%2
                  y0 = (1 + dimy)%/%2
                  z0 = (1 + dimz)%/%2
                  x1 = x0
                  y1 = y0
                  z1 = (1 + dimz)%/%4
                  x2 = x0 - (1 + dimx)%/%4
                  y2 = y0
                  z2 = z0
                  x3 = x0 + r_motor
                  y3 = y0 + r_motor
                  z3 = z0 + r_motor
                  
                  if ((x - x1)^2 + (y - y1)^2 + (z - z1)^2 <= r_motor^2 | (x - x2)^2 + (y - y2)^2 + (z - z2)^2 <= r_motor^2 | (x - 
                    x3)^2 + (y - y3)^2 + (z - z3)^2 <= r_motor^2) {
                    motor_3d[x, y, z] = 0
                  }
                }
            }
        }
        
        # mask is not provided
    } else {
        
        r_mask = (min(dimx, dimy, dimz) - 2)%/%2
        mask = array(0, dim = dim_data)
        for (x in 1:dimx) {
            for (y in 1:dimy) {
                for (z in 1:dimz) {
                  x0 = (1 + dimx)%/%2
                  y0 = (1 + dimy)%/%2
                  z0 = (1 + dimz)%/%2
                  x1 = x0
                  y1 = y0
                  z1 = z0 - (1 + dimz)%/%4
                  x2 = x0 - (1 + dimx)%/%4
                  y2 = y0
                  z2 = z0
                  x3 = x0 + r_motor
                  y3 = y0 + r_motor
                  z3 = z0 + r_motor
                  if ((x - x0)^2 + (y - y0)^2 + (z - z0)^2 <= r_mask^2) {
                    mask[x, y, z] = 1
                  }
                  if ((x - x1)^2 + (y - y1)^2 + (z - z1)^2 <= r_motor^2 | (x - x2)^2 + (y - y2)^2 + (z - z2)^2 <= r_motor^2 | (x - 
                    x3)^2 + (y - y3)^2 + (z - z3)^2 <= r_motor^2) {
                    motor_3d[x, y, z] = 0
                  }
                }
            }
        }
        
    }
    
    tspan = ons[length(ons)] + 2 * dur[length(dur)] - 1
    data_4d = runif(dimx * dimy * dimz * tspan, 0, 15)
    dim(data_4d) = c(dimx, dimy, dimz, tspan)
    t_range = c(1:tspan)
    # on_time = t_range[rep(c(TRUE, FALSE), c(10, 10))] off_time = t_range[!rep(c(TRUE, FALSE), c(10, 10))]
    
    
    on_time = c()
    off_time = c()
    # length of ons and dur should have same length
    for (i in 1:length(dur)) {
        
        on_st = ons[i]
        on_end = ons[i] + dur[i] - 1
        on_time_pts = c(on_st:on_end)
        off_st = ons[i] + dur[i]
        off_end = ons[i] + 2 * dur[i] - 1
        off_time_pts = c(off_st:off_end)
        on_time = c(on_time, on_time_pts)
        off_time = c(off_time, off_time_pts)
    }
    
    for (t in t_range) {
        data_4d[, , , t] = data_4d[, , , t] * mask
        if (t %in% on_time) {
            on_motor_data = runif(dimx * dimy * dimz, 12, 18)
            dim(on_motor_data) = dim_data
            data_4d[, , , t] = data_4d[, , , t] * motor_3d + on_motor_data * (1 - motor_3d)
        }
    }
    return(list(fmri_data = data_4d, mask = mask, ons = ons, dur = dur, on_time = on_time))
}